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liYTRODUCTION 

A  time  series  can  be  defined  as  a  sequence  of  discrete  observations 
taken  at  uniform  intervals  on  a  time  scale.   Since  time  series  occur 
frequently  in  industrial  and  economic  situations  and  since  many  decisions 
are  based  on  the  predictions  of  the  future  values  of  these  time  series  it 
is  an  important  decision  aid  to  develop  a  method  for  forecasting  these  time 
series. 

In  order  to  treat  the  time  series  analytically  it  is  described  as  the 
sum  of  two  components.   The  first  is  the  process  which  generates  the  series 
and  the  second  is  some  superimposed  random  noise,  or  variation  in  the  obser- 
vations.  Hence,  the  time  series  can  be  represented  as 

y(t)  ■  i(%)   +  e(t) 

where  £(t)  is  the  value  of  the  process  at  time  t  and  e(t)  is  the  random 
noise  associated  with  the  observation  at  time  t.   The  distribution  of  the 
noise  samples  has  the  following  properties: 

(1)  The  expected  value  is  zero,  i.e. 

E(e(t))  =  0 

(2)  The  noise  samples  have  no  series  correlation,  i.e. 

E(e(i)e(j))  =0  for  i  4   j 
Considering  the  nature  of  the  time  series  the  object  of  the  forecasting 
technique  is  to  provide  a  forecast  that  will 

(1)  Dampen,  or  smooth,  the  effect  of  the  random  noise, 

(2)  Reflect  any  trends  which  the  time  series  might  be  undergoing, 
i.e.  the  time  series  need  not  be  stationary, 

and 

(3)  Provide  an  unbiased  estimate  of  the  time  series. 


Hence,  a  forecasting  technique  is  sought  which  will  seek  a  balance  between 
the  response  to  secular  trends  in  the  process  and  the  errors  caused  by  the 
random  noise. 

The  process  which  generates  the  time  series  can  be  described  by  two 
components.   These  are  the  secular  trend  of  the  mean  and  the  periodic  com- 
ponent of  the  series.   Along  with  the  trend  which  the  mean  of  the  time 
series  might  be  undergoing  the  periodic  component  might  also  contain  trends. 
These  latter  trends  are  of  two  types 

(1)  Changing  amplitudes 

(2)  Shifting  phase  angles 

Therefore,  the  forecasting  model  must  adequately  represent  the  time  series 
and  be  able  to  adjust  to  these  trends.   It  is  found  that  forecasting  models 
conposed  of  polynomials ,  trigonometric  functions  and  multiples  of  these 
functions  fulfill  the  requirements  of  the  forecasting  model. 

These  forecasting  models  can  be  represented  in  vector  form  as 
y(t)  =  a'f(t) 
where  a  is  a  vector  of  coefficients  and  f(t)  the  vector  of  fitting  functions 
evaluated  at  time  t.  The  requirements  of  the  forecast  are  satisfied  in  the 
method  by  which  the  vector  of  coefficients  is  estimated.   The  estimation  of 
this  vector  is  performed  by  the  use  of  general  exponential  smoothing.   In 
this  process  the  estimate  of  the  coefficients  is  based  on  the  discounted 
least  squares  criterion.   That  is,  the  sum 

I   BJ(y(-t)  -  y(-t-l)) 

is  minimized  where  y(t)  is  the  value  of  the  time  series  at  time  t. 


y(t-l)  is  the  forecast  of  this  value  at  the  previous 
period, 
8  is  a  positive  constant  less  than«  one . 
This  forecasting  technique  has  the  following  properties: 

(1)  The  weight  given  the  observations  in  the  forecast  is  discounted 
on  a  time  scale  at  a  rate  which  can  be  controlled  by  the  value 
of  the  discount  factor  B,  thereby  yielding  control  of  the  re- 
sponsiveness of  the  forecast. 

(2)  All  past  data  is  contained  in  a  single  word  of  information. 

(3)  A  simple  recursive  relationship  can  be  developed  for  re-evaluating 
the  estimates  of  the  vector  of  coefficients  with  each  observation. 

The  justification  for  the  use  of  general  exponential  smoothing  lies 
in  the  proof  that  if  the  forecasting  model  is  a  true  representation  of  the 
process  then  the  expected  value  of  the  vector  of  coefficients  obtained  by 
exponential  smoothing  will  be  the  true  value  of  this  vector.  The  problem, 
then,  becomes  that  of  determining  the  forecasting  model.   It  is  found  that 
the  trend  of  the  mean  can  be  represented  by  a  polynomial.   To  this  end 
polynomial  regression  is  used  to  fit  the  curve  representing  the  trend  of 
the  mean  to  the  data.   Having  performed  this  regression  both  qualitative 
and  quantitative  estimates  for  the  terms  in  the  forecasting  model  which 
represent  the  trend  of  the  mean  are  available.   The  regression  curve  is  then 
subtracted  from  the  data,  an  operation  referred  to  as   detreuding  • 

Making  the  data  available  in  the  "detrended"  form  is  the  first  phase  in 
the  analysis  of  the  periodic  component.   I&aentially  "detrending"  allows 
the  periodic  component  to  be  observed  sepernlely .   The  second  analysis  per- 
formed on  the  periodic  component  is  autocorrelation  analysis.   The  auto- 
correlation function  is  defined  and  shown  to  have  a  maximum  at  the  basic 


period  of  the  data.   Moreover,  autocorrelation  analysis  can  be  used  as  a 
test  of  significance  that  a  process  exists. 

The  basic  theorem  of  Fourier  series  states  that  if  a  discrete  series 
is  basically  periodic  then  the  coefficients  of  a  series  of  sines  and  cosines 
of  the  fundamental  harmonics  of  that  period  can  be  determined  so  that  the 
Fourier  series  yields  the  data  points.  This  analysis  is  extended  to  a 
trigonometric  series  which  includes  all  the  frequencies  in  the  data  and  then 
to  a  least  squares  evaluation  of  these  coefficients.   It  is  shown  that  a 
measure  of  the  contribution  of  each  frequency  in  describing  the  time  series 
can  be  expressed  as  a  function  of  these  coefficients  and  that  using  this 
measure  the  periodic  component  can  be  adequately  expressed  as  a  limited 
number  of  these  frequencies. 

A  study  is  made  of  the  sensitivity  of  the  forecasts  obtained  through 
the  use  of  exponential  smoothing.  The  sensitivity  analysis  is  carried  out 
on  the  following  parameters 

(1)  The  fitting  functions  used  to  describe  the  forecasting  model, 

(2)  The  basic  period  of  the  forecasting  model,  and 
13)  The  value  of  the  discount  factor. 

It  is  found  that  the  choice  of  each  of  these  parameters  significantly  affects 
the  accuracy  of  the  forecasts.   In  the  case  of  the  smoothing  constant  a 
literature  search  reveals  that  no  method  of  determining  the  optimal  value  of 
this  constant  has  been  presented.  This  research  does  not  attempt  to  develop 
such  a  method.   It  is  felt,  however,  that  a  local  optimal  value  for  the 
discount  factor  might  exist  for  a  particular  time  series  and  set  of  fitting 
functions.   Hence,  a  parametric  investigation  of  the  smoothing  constant  is 
carried  out. 


A  general  program  to  perform  exponential  smoothing  was  written.   This 
program  has  the  ability  to 

(1)  Change  the  forecasting  model, 

(2)  Change  the  basic  period  of  the  forecasting  model, 

(3)  Change  the  value  of  the  smoothing  constant  and 
(U)   Change  the  time  series. 

This  program  has  the  dual  purpose  of  serving  as  a  medium  for  carrying  out 
and  evaluating  the  effectiveness  of  general  exponential  smoothing  for  fore- 
casting a  time  series  and  determining  the  feasibility  and  economy  of  per- 
forming this  forecasting  technique  using  a  Fortran  type  of  processor. 
Moreover,  the  results  can  be  compared  with  those  obtained  by  Brown  (l) 
to  verify  the  consistancy  and  accuracy  of  the  program.   The  general  nature 
of  the  program  is  necessary  in  order  to  carry  out  the  sensitivity  analysis 
on  the  forecast  parameters  and  the  parametric  investigation  of  the  smoothing 
constant.   Furthermore  this  type  of  general  program  can  serve  to  illustrate 
that  the  requirements  of  an  industrial  situation  wherein  many  and  varied 
time  series  are  encountered  can  be  satisfied. 

Using  the  I.B.M.  lUlO  system  at  Kansas  State  University  the  program 
for  general  exponential  smoothing  can  not  be  accomodated  in  the  available 
core  storage  capacity.   Hence,  phasing  of  the  program  is  necessary.   The 
system  at  Kansas  State  University  is  programmed  internally  with  PR-155  and 
has  seven  magnetic  tape  drive  units  which  makes  phasing  possible.   Each 
phase  is  run  independently  and  upon  its  completion  the  processor  auto- 
matically clears  core  and  loads  in  the  next  phase.   Any  information  required 
for  following  phases  is  written  onto  a  work  tape.   In  this  system  one  work 
tape  is  available,  however,  after  the  program  is  compiled  two  more  of  the 


tapes  can  be  used  as  scratch  files.  These  three  tape  units  are  the  min- 
imum required  to  perform  the  internal  data  transmission  between  phases. 

Two  time  series  are  used  in  the  application  of  the  forecasting 
techniques.  The  first  is  the  number  of  miles  traveled  on  international 
airline  routes  measured  at  monthly  intervals  from  January  1949  to  December 
I960.   This  data  was  obtained  from  Smoothing  Forecasting  and  Prediction  of 
Discrete  Time  Series  by  Robert  Goodell  Brown.   The  second  time  series  used 
was  the  sheep  population  in  England  and  Wales  measured  in  yearly  intervals 
from  1867  to  1939-  This  data  was  obtained  from  The  Advanced  Theory  of 
Statistics  Volume  II  by  M.  G.  Kendall. 


1.   REPRESENTATION  OF  THE  TIME  SERIES 

1.1  The  Time  Series  Model 

A  time  series  is  a  sequence  of  observations  taken  at  equal  time 
intervals.  For  the  purpose  of  analysis  the  time  series  can  be  considered 
to  be  made  up  of  two  elements 

(1)  The  process  which  generates  the  time  series 

(2)  Some  superimposed  random  noise 

Thus  the  time  series  may  be  represented  on  the  following  manner 
y(t)  =  £(t)  +  e(t) 

where  5(t)  is  the  process 

e(t)  is  the  noise  in  the  t   observation. 

The  distribution  of  e(t)  has  the  properties 
EUt)  =  0 

E(ete  )  =0  for  i  ?   J 
=  o2  for  i  =  J 

where  o2  is  the  variance  of  the  noise  distribution. 
Recognizing  that  the  noise  in  the  time  series  is  random,  no  attempt  is  made 
to  represent  it.  The  techniques  of  representing  the  time  series  are  con- 
cerned with  the  process.  Due  to  the  conditions  which  generate  time  series 
it  is  generally  adequate  to  describe  the  process  in  terms  of  two  components 

(1)  The  trend  which  the  mean  of  the  series  is  following, 

(2)  A  cyclical  component  which  is  superimposed  upon  this  trend. 
The  trend  component  may  be  represented  by  the  following  functions 

(1)  Polynomials, 

(2)  Exponentials, 
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whereas  the  periodic  component  must  necessarily  be  represented  by 

trigonometric  functions.  This  representation  is  based  on  the  Fourier 

analysis  which  deals  with  functions,  either  continuous  or  discrete,  by  means 

of  a  series  of  fundamental  harmonics.  The  principal  theorem  of  Fourier  series 

may  be  stated  as  follows: 

If  f(t)  is  a  single-valued  function  which  has  a 
derivative  throughout  the  interval  -a<t<a  except  for 
a  finite  number  of  points  at  which  it  has  finite  dis- 
continuities, and  for  other  values  of  t  is  defined  by 
the  equation 

f(t  +  2a)  =  f(t) 

then  f(t)   can  be  represented  by  means  of  the  Fourier 
series 

y  =  !jA     +  A     cos(nt/a)   +  A     cos(2Ht/a)   + 

+  A     cos(3Ht/a)   +    . . . 
+  B.    sin(Ht/a)   +  B     sin(2JIt/a) 

+  B     sin(3nt/a)  +   ...  (1:1:1) 

See  appendix  A  for  the   development   of  the  Fourier  analysis. 

In  the  application  of  the  Fourier  series  to  representing  the  periodic 
part  of  the  time  series   only  the  terms   in   (1:1:1)  which  are  shown  to  be 
significant  will  be  used.     A  limited  number  of  terms   are  necessary  to  strike 
a  balance  between  the  accuracy  of  the  model  and  the  length  of  the   compu- 
tations.     Hence  a  general  representation  of  the  time   series   is: 

y  =  trend  +  A(T)    cos(2nt/T)    +  B(T)    sin(2Ut/T) 

(1:1:2) 
+  A(T' )    cos(2nt/T' )   +  B(T')    sin(2nt/T') 

where  T  and  T'  are  the  periods  of  the  harmonics  which  show  significant  con- 
tribution to  the  representation  of  the  time  series. 


The  general  model,  (1:1:2),  is  adequate  for  representing  a  time  series 
which  displays  a  secular  trend  in  the  mean  and  a  periodic  component  super- 
imposed on  that  trend.  In  some  time  series,  however,  two  other  types  of 
trends  may  appear  in  the  periodic  component 

(1)  Shifting  phase  angles, 

(2)  Growing  amplitudes. 

The  technique  for  representing  these  trends  is  to  include  in  the  model  a 
set  of  terms  of  the  form 


(a  +  at)  cos(2Ht/T)  +  (a  +  a^t)  sin(2IIt/T) 


(1:1:3) 


where  T  is  the  harmonic  in  which  either  or  both  of  the  above  trends  occur. 
An  example  of  the  representation  of  a  time  series  might  be  sighted  as 
follows.  Consider  the  time  series  in  figure  (1:1:1) 


1949     1950     1951     1952     1953     1954     1955     1956     1957     1956     1959     i960 


Fig.  1:1:1  International  Airline  Passenger  Data 

Looking  ahead  in  the  analysis  momentarily  this  time  series  might  be  ade- 
quately represented  by  a  linear  trend  with  a  periodic  component  of  12  months. 
Moreover,  and  without  any  Justification  at  this  point,  a  harmonic  at  o  months 
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might  be  significant.  Furthermore,  the  amplitude  of  the  periodic  com- 
ponent seems  to  be  increasing.  With  these  ideas  the  terms  included  in 
the  model  would  he: 

(1)  a  +  at     to  represent  the  linear  trend, 

(2)  (a  +  a^t)  sin(2nt/12)  +  U3  +  a?t)  cos(2nt/12) 

to  represent  the  12  month  periodic  component  with 
increasing  amplitude, 

(3)  a6  sinUnt/12)  +  a?  cos(l*Ht/12) 

to  represent  the  harmonic. 

1:2  Trend  Analysis  of  a  Time  SerieB 

The  previous  section  showed  that  an  adequate  representation  of  the 
time  series  depended  on  proper  evaluation  of  the  trends  that  the  series 
was  following.  These  trends  are  of  two  types: 

(1)  Trend  of  the  mean 

(2)  Trend  of  the  periodic  component. 

Moreover,  the  analysis  of  the  periodic  component  was  treated  independently 
of  the  time  series.  In  effect  this  independent  treatment  is  equivalent  to 
evaluating  the  series  with  the  trend  of  the  mean  removed.  The  purpose  of 
trend  analysis  of  the  time  series,  then,  is  twofold. 

(1)  Trend  analysis  provides  a  quantitative  estimate  of  the  trends 
in  the  time  series . 

(2)  Trend  analysis  removes  the  trend  of  the  mean  from  the  time 
series.   This  "detrended"  form  of  the  data  is  used  in  further 
analysis . 
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The  mechanism  used  for  trend  analysis  is  polynomial  regression. 

(see  appendix  B).  In  the  expression  for  time  time  series 

y  =  6(t)  +  e(t) 

the  process  S(t)  can  be  written  as 

£(t)  =  trend  +  periodic  component 

If  the  trend  component  is  expressed  as  y"  then  the  expression  for  the  time 

series  with  the  trend  removed  is 

y'  =  C(t)  +  e(t)  -  y" 

=y_y"  (1.2.1) 

where  y"  is  the  value  of  the  regression  curve.  Hence  y'  is  the  periodic 

component  of  the  time  series. 

In  selecting  the  regression  curve  to  represent  the  trend  of  the  mean, 

care  must  be  taken  not  to  include  any  of  the  periodic  component  in  the  trend 

to  be  removed.  The  trigonometric  functions  can  be  expressed  as  a  series 

3         5         7 
tJ   .    tJ       V 


3int  =  t-2T+5T-yr+    ••• 


t2       tU       t6 

cost  =  i-2T+-^r-^r  + 


(1.2.2) 


Now  suppose  the  regression   curve   is   taken  to  be  a  polynomial  of 
degree  k. 


2  k 

y   =   a     +   a  t  +   a  t      +    . . .    +   &t  (1.2.3) 
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If,   in  fact,   the  time  series   consists  of  a  trend  component  which  can  be 
expressed  as   a  polynomial  of  degree  k',  where  k'    <  k,    and  a  periodic  com- 
ponent which  can  be  expressed  as 

y'   =  A(T)   cos(2nt/T)   +  B(T)   sin(2JIt/T) 

then  considering   (1:2:2)   the  series     y     can  be  expressed  as  a  polynomial 
of  degree  k" ,   where  k"    >  k. 

y  =  aQ  +  &1t  +  a2t2  ♦....♦  a^t*" 

Hence  if  the  order  of  the  regression  takes  on  the  value  k,  where  k  >  k' 
then  the  regression  curve  obtained  can  be  expressed  as  the  sum  of  two  series 

ya  =  a0  +  "lt  +  a2t2+    •••   +  \*k 

y8  =   60+   V+82t2+    "•   +   Vk 

where  y     is  the  series  contributed  by  the  secular  trend  of  the  mean  and  y 

6 

is  the  series  contributed  by  the  periodic  component.  Hence,  the  value 
obtained  for  the  regression  curve  will  be 

y"  =  y"  +  y" 
where  in  fact  it  should  be 


Therefore  if  k  >  k'  then  part  of  the  periodic  component  will  be  removed 
from  the  data  in  "detrending"  and  the  trend  analysis  will  give  a  distorted 
picture  of  the  data.   That  is  the  "detrended"  data  will  be 

y'  =  y  -  (y0  ♦  yg) 

instead  of 
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r  =  y  -  ya  ■ 

An  example  of  trend  analysis  can  be  found  in  its  application  to  the 
time  series  in  figure  (1:1:1).  In  this  case  the  secular  trend  of  the 
mean  is  taken  to  be  a  linear  relationship.  That  is  the  expression  for  the 
trend  of  the  mean  is  taken  to  be 

y"  =  a0  +  V 

where  a     and  a     are  the  coefficients   to  be  determined  by  regression  analysis. 
Appendix  B  gives  the  expressions   for  the  least  squares  estimate  of  these 
coefficients  as 


n  n 

I     ty  -  I     ly 

i=l 1=1 


n       „  n  n         „ 

I     t2  -  2t      I  t  +      I    (t)2 


1-1  i-1  i-1 

Hence  the  expression  for  the  detrended  time  series  is 

y'(t)  =  y(t)  -  t  Bj 

where  t  denotes  the  t   value  in  the  time  series 
Thus  far  the  trend  analysis  has  fulfilled  the  purposes  of  providing 
the  time  series  in  detrended  form  and  providing  qualitative  estimates  for 
the  secular  trend  in  the  mean.  A  third  purpose  of  trend  analysis  is  that 
of  providing  quantitative  estimates  for  the  trends  of  the  periodic  component 
of  the  time  series.   With  the  time  series  expressed  in  the  form 


lit 


y'  =  y  -  y" 

the  periodic  component  can  be  looked  at  directly  without  any  effect  of  the 
secular  trend  of  the  mean.  Thus  if  the  amplitude  of  the  periodic  component 
is  undergoing  a  secular  trend  the  series  would  appear  in  the  detrended  form 
as 


,-ff^o" 


°oo  o° 


Fig.  1:2:1  Example  of  "detrended"  time  series 

in  which  case  the  slope  of  the  line  A'A'  can  be  easily  determined.  If 
however,  the  mean  is  following  a  secular  trend  which  could  be  expressed  as 


2 

0  '  "l"  '  ~2" 


y"  =  a„  +  a,  t  +  a„t 


the  series  before  trend  analysis  is  applied  might  look  like  figure  (1:2:2). 
Hence,  the  preceeding  analysis  would  be  much  more  difficult  to  perform. 
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o   o 


o    o 
o 


°°OoO 


°o  oo 


Fig.  1:2:2  Series  Before  "Detrending" 

1.3  Autocorrelation  Analysis  of  the  Time  Series 

In  the  previous  section  the  time  series  was  expressed  as 
y(t)  =  C(t)  +  e(t) 

where  £(t)  is  the  process 

e(t)  noise  samples,  where  e(t)  is  the  t   sample 

from  a  probability  distribution  with  zero  mean 
The  application  of  the  principal  theorem  of  Fourier  series  is  based 
on  the  assumption  that  a  periodic  component  of  the  time  series  exists  which 
can  be  represented  as 

f(t  +  2a)  =  f(t) 

The  purpose  of  autocorrelation  analysis  is  to  answer  three  questions 
about  the  time  series 

(1)  Does  a  process  exist? 

(2)  Does  a  periodic  component  exist? 

(3)  What  is  the  basic  period  of  the  time  series? 
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The  mechanism  used  in  answering  these  questions ,  correlation  analysis ,  is 
developed  in  Appendix  C. 
Consider  the  term 

T 

I     r.y.k  (1.3.1) 

j=k+i  J  J  K 

where  y  is  the  J   observation  in  the  time  series 

j  =  1,2,...,  T 

k  is  the  lag  between  the  terms  in  the  summation 
The  expected  value  of  this  term  is  denoted  as  the  average  lagged  product. 
In  particular,  if  the  sequence  y  has  been  adjusted  so  that  the  expected 
value,  E(y)  =  0,  then  the  average  lagged  product  is 
T 

the  autocovariance. 

The  variance  of  a  sequence  of  numbers  that  have  been  adjusted  to  have 

a  mean  value  of  zero  i3  Just  the  expected  value  of  the  squares  of  these 

o 
numbers.   This  can  be  expressed  as  E(y  )  which  is  effectively  R   (0). 

Moreover,  the  autocovariance  expressed  in  normalized  form 

p(k)  =  Rxx(k)/Rxx(0)  (1.3.10 

is  merely  the  autocorrelation  coefficient.   The  set  of  values  for  the  auto- 
correlation coefficient  for  all  lags, 
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k-  1,  2,  

is  defined  as  the  autocorrelation  function.  Hence,  the  autocorrelation 
function  for  a  lag  k=0  is  one.  Three  other  properties  of  the  autocor- 
relation function  are  of  significance  in  the  analysis. 

(1)  The  range  of  p  is  between  +  1. 

(2)  Pure  random  noise  will  have  zero  correlation  between  samples 
not  identically  equal  to  each  other. 

(3)  If  whenever  y   is  positive  so  is  y   .  And  whenever  y^  is 
negative  so  is  yt+v-  Then  the  autocovariance  will  be  large 
and  positive.  In  this  case,  pairs  of  observations  k  units 

of  time  apart  in  the  sequence  are  highly  correlated  and  one  can 
be  used  to  forecast  the  other.  In  a  similiar  manner  if  y 
positive  usually  implies  that  yt+k  is  negative  (and  vice  versa), 
one  can  still  be  used  to  forecast  the  other.  The  autocovariance 
in  this  case  will  be  large  and  negative  (l,391*). 
In  the  application  of  autocorrelation  analysis  Brown  (1,395)  suggests 
the  following  proceedure.   Plot  the  observations  to  see  if  you  should 
expect  a  secular  trend  or  a  significant  cyclical  pattern.  If  there  is  a 
secular  trend,  fit  a  straight  line  to  the  data  by  least  squares.  Using  this 
least  squares  fit  adjust  the  data  to  zero  expected  value.   Compute  the 
autocovariances  using  (1:3:2). 

Rxx(k)  =  (l/T-l-k)  jk+iVj_k 

This  suggested  method  of  Browns  can  be  applied  more  discretly  by 
taking  advantage  of  the  trend  analysis  in  the  previous  section.   In  that 
section  the  data  was  obtained  in  the  detrended  form  as 
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y    =  y  -  y 

where     y     is   the  time   series 

y"   is  the   secular  trend  of  the  mean 
Moreover,    since   the  time  series   can  be  expressed  as 

y  =  trend  +  periodic   component  +  random  noise 
it   follows   that 

y'    =  periodic   component  +  random  noise 
It  was   shown  in  appendix  A  that  by  the  principal  theorem  of  Fourier 
series   the  periodic   component   can  be   represented  by  the   series 

y1    =  hk     +  A     cos(nt/a)   +  A     cos(2nt/a)   +  A     cos(3nt/a)   +    ... 

+  B,    sin(nt/a)    +  B     sin(2IIt/a)   +  B     sin(3nt/a)   +    ... 
Note  that  the  expected  value  of  this   series   is  merely  h&~-     If  this 
expected  value  is  subtracted  from  the   "detrended"   series  the   following  is 
obtained. 

y  =  y'   -  hkQ  +  e(t) 

where  e(t)  is  the  noise  sample 
The  expected  value  of  this  series  can  be  expressed  as 

K(y)  =  E(y'  -  isAj  +  E(e(t)) 
Now  since  the  expected  value  of  the  first  term  on  the  right  was  shown  to 
be  zero  and  since  by  definition  the  expected  value  of  the  noise  is  zero, 
then 

K(y(t))  =  0 
and  this  series  satisfies  the  conditions  for  the  application  of  (1.3.2), 
the  autocovariance . 
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'  To  facilitate  the  use  of  the  autocorrelation  function  in  answering 
the  question  previously  posed  in  this  section,  i.e. 

(1)  Does  a  process  exist? 

(2)  Does  a  periodic  component  exist? 

(3)  What  is  the  basic  period  of  the  series? 
some  of  the  properties  of  this  function  must  be  noted. 

An  example  of  autocorrelation  analysis  is  given  bt  H.  T.  Davis 

+  1.00 

+  0.50 

0 

-0.50 

- 

+  0.50 

0 

-0.50 
-  1.00 

- 

/ 

- 

"'•1?70-r60 ""~3o      -»>   "    0      W0      »»      ♦«!  ♦'" 

Figure  1:3:1  Autocorrelation  function  for  industrial  stock 
prices,  t  measured  in  months 

Referring  to  this  plot  Davis  comments: 

'  "It  will  be  observed  from  the  graph  that  the  function 
damos  rapidly.   It  changes  from  positive  to  negative  at  ap- 
proximately t  =  +  10,  and  again  becomes  positive  at  t  -  +  W. 
As  we  shall  show-later  on  this  may  be  interpreted  as  indicating 
a  cycle  of  1*0  months."  (3.356) 

Hence  the  autocorrelation  function  will  reach  its  maximum  at  the  basic 

period  of  the  series.  The  above  characteristic  may  be  proved  as  follows. 
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The  autocovariance  (1.3.2)  may  tie  expressed  as 

T 

,k)  =  11m  — — - 
xx 


1 
Rxx(k)  =  llm  2T  *  1  t   y(t+k)  y(t) 


t-*»       t=-T 

where  the  series  y(t),  t  =  1,2,...  T,  has  been  adjusted  so  that  the  expected 
value  is  zero.   But  the  series  adjusted  in  this  manner  can  be  expressed  as 

y'(t)  =  y(t)  -  y"(t) 
where  y"(t)  is  the  secular  trend  of  the  mean.   Therefore,  the  adjusted 
series  can  be  expressed  as 

y'(t)  =  periodic  component  +  noise 
In  Appendix  A  it  is  shown  that  the  periodic  component  can  be  expressed  as 
a  Fourier  series 

n  n 

y'(t)  =  a„  +  )   a.  cos  id.  +  ;  b.  cos  id. 

0   .  <•,   i      i   .  L,      i      i 

1=1  i=l 

Making  the  transformation  from  the  general  formulation  of  the  Fourier 
series  above  to  an  infinite  series  of  cosines  (this  change  will  only 
simplify  the  calculations),  the  autocovariance  may  be  re-expressed  as 

,    T    n   n 
Rxx(k)  ■  £lm2TTI  j_T  l=1     ^Vj  "■(-!*]  «•(•!<«♦*>]  +  Vt+k 

since  the  cross  product  between  the  noise  and  the  cosine  signal  have 
expected  value  zero.  The  expected  value  of  all  terms  of  the  form 

cos(w  t)  cos  (id  t) 

is    also   zero   for  i    4   J.      Therefore,    ttie   autocovariance   reduces    to 


21 


R  (k)  ■  h     7  a.  cos  w.k  +  R  (k) 
xx        *  1      1     ee 


If  the  assumption  is  made  that  the  noise  has  no  serial  correlation  then 

REE(k)  =  o*6(k) 

R     (k) 
and  the  autocorrelation  function  p(k)   =  j—nrr ;     will  have   a  local  maximum 

xx 

at  R  =  (2Jl/u.  ).   This  completes  the  proof.  (1,396) 

The  autocorrelation  function  can  be  used  to  determine  if  a  process 
exists.   Davis  (Jjl'O)  provides  an  example  of  the  autocorrelation  function 
for  a  completely  random  series.   The  random  series  was  constructed  in  the 
following  manner.  The  percentages  of  trend  of  the  Dow-Jones  industrial 
averages  for  the  prewar  period  (1897-1913)  were  written  on  cards  and  these 
cards  were  drawn  at  random  to  form  a  series  of  20i»  items,  that  is  N  '=  20I4. 
The  polt  of  the  autocorrelation  function  was  determined  to  be: 


Fig.  1:3:2  Autocorrelation  of  a  random  series-  the  dotted 
lines  define  the  standard  error  band 
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The  standard  error  band  is   computed  by  the  relationship 

S2     .  gi  (82  _  *2  s2) 

yx       U-2       y  x 

where  b  is  given  by  the  regression  formula  (5) 

b  =  . 

x2  -  ([  x.)2/N 

The  standard  error  band  varies  from  +0.070  at  the  beginning  vhere  N  =  204, 
to  +.0076  at  the  end,  where  H  =  20k  -  30  =  Ilk.     Hence  the  distribution 
of  the  lagged  values  is  seen  to  meet  the  test  of  randomness  in  a  satis- 
factory manner. 

1:U  Spectral  Analysis  of  the  Time  Series 

The  original  definition  of  the  time  series  took  the  form 

y  =  t(t)  +  e(t) 

where  t(t)  is  the  process 

e(t)  is  the  random  noise 
Moreover,  the  process  £(t)  can  be  expressed  as 

t(t)  =  trend  +  periodic  component 
The  purpose  of  this  section  is  to  investigate  the  periodic  component  of 
the  series . 

The  principal  theorem  of  Fourier  series  (appendix  A)  shows  that  the 
periodic  component  of  the  time  series  can  be  represented  to  any  desired 
accuracy  by  the  series 
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y'   =  ijA0  +  A±  cos(nt/a)  +  A2  cos(2nt/a)   +  A     cos(3nt/a)   +   ... 
+  B1  sin(nt/a)   +  Bg  sln(2nt/a)   +  B     sin(3nt/a)   +   ... 

(l:lt:l) 
The  terms  of  the  above  series   represent  the  harmonics   of  the  basic 
period  of  the  time  series.      In  the  general  case  the  process   can  be  represented 
more  accurately  by  increasing  the  number  of  terms   in  this   series.     However, 
these  harmonic  terms  together  with  the  terms  which  represent  the  secular 
trend  of  the  mean  form  the  forecasting  model.     Hence,   as  the  number  of  terms 
in  the  Fourier  series  increases  the  time  series   is   discribed  more  accurately 
but  the   calculations  for  the  forecast  are  also  increased. 

The  study  of  harmonic  analysis   shows  that  the   frequencies   represented 
in   (l:U:l)   differ  in  their  contribution  to  representing  the  time  series. 
The  purpose  of  the  spectral  analysis  is  to  obtain  a  measure  of  the  contri- 
bution of  each  frequency.     This  measure  is  used  as   a  basis   for  selecting 
the  frequencies  of  the  periodic  terms   to  include  in  the   forecasting  model. 

The  analysis  of  the  representation  of  the  periodic  component  of  the 
time  series  by  the  Fourier  series   (lsbtl)   is  based  on  the  variance  of  the 
approximation.      In  this   analysis  the  inequality  of  Bessel  is   used  to  express 
the  variance  of  (l:lt:l)   in  terms  of  the  Fourier  coefficients. 

In  order  to  derive  this  inequality  assume  that  the  process  y'(t)  has 
been  approximated  by  the  first  N  harmonics  of  the  Fourier  series  (1:U:1), 
that  is 

N  N 

y'(t)  *  hkQ  +      [  A     cos(nnt/a)   +      [  B     sln(nllt/a)  (l;l»-2) 

n=l  n=l  n 


2k 


If  the  right-hand  member  of   (l:U:2)   is  represented  by  yn(t)   and  the 
intergal  of  the  square  of  the  residual  is   considered,   then 

2 

I  =  i/a    (y'(t)   -  y   (t))      dt 
a  -a    w  n       ' 

By  expansion 

=  a"  £    ^y'2(t)   _  *,<*^n(t)   +  ynU)^    " 
Taking  account  of  the  integrals 


/a  sin(mnt/a)   sin(nllt/a)   dt  =  /     cos(mnt/a)   cos(nHt/a)   =  0 


for  m  4  n 


-  !&  sin2(nnt/a)   dt  =  -  f&  cos2(nnt/a)   dt  =  1 
a  -a  a  -a 


/a  sin(mnt/a)   cos(nHt/a)   dt  =  0, 
-a 

and  observing  the  definitions  of  the  Fourier  coefficients  given  in 
Appendix  A.     The  expression  for  the  integral  I   can  be  obtained  as 

*  =  zil*'*M  dt  -  {hAl +  Ri +  R2 +  ■  ■  -  +  4] 

2  2  2 

where  R     =  A     +  B 
n  n  n 

Moreover,  since  the  integrand  of  the  integral  is  positive  or  zero  the  integral 

itself  is  positive  or  zero,  and  thus  the  Bessel  inequality  for  Fourier 

coefficients  is  obtained.  (3,65) 
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^0  +  Rl  +  R2+  •  •  •  +  4   "  ha(f(t»2  dt  (1«*«SJ 

By  noting  that  the  arithmetic  average  of  y'(t)  is  equal  to  %k   ,  the 


At  =  y'(t)  -  yn(t) 
is  given  by 

O   -  *C  -  H     I    (Hf)  -  h     \(t?  +  J?)  (a:U:U) 

n=l        n=l 

The  term  which  is  used  as  a  measure  of  the  contribution  of  each  harmonic 
is 

2o2 

where  R2(T)  =  A2(T)  +  B2(T) 

T  is  the  period  of  the  harmonic 

a2   is  the  variance  of  the  data. 
From  Bessel's  theorem,  the  variance  a2,   of  the  series  after  n  terms 
have  been  removed;  that  is,  equivalent^,  if  the  series  were  corrected  for 
these  harmonics,  is  (3,71) 

0*  -  (1  -  [En)02 

Hence,  the  energies  of  the  harmonics  as  expressed  by  (l:U:5)  are  strictly  ad- 
ditive if  the  harmonics  belong  to  the  Fourier  sequence.  If  the  harmonics  do 
not  belong  to  the  Fourier  sequence  then  this  expression  is  only  approximately 
correct.  Hence,  E(T)  is  an  appropriate  measure  of  the  contribution  of  a 
harmonic. 
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The  following  expressions  for  the  Fourier  coefficients  for  the  case  of 
a  discrete  series  are  given  in  Appendix  A. 

?  T 
A  =  -  I   f  cos(2nnt/T) 

n   X  t=l  * 

(1.1.. 6) 

?  T 
B  =|  I   f  sin(2nnt/T) 

"     t»l 

where  T  is  the  number  of  observations  in  the  period  of 
the  harmonic 
f  is  the  time  series  evaluated  at  time  t. 
Upon  consideration  of  these  formulas  a  major  objection  is  noted.   The 
frequencies  evaluated  by  equations  (l.U.6)  reflect  only  those  frequencies 
which  are  present  in  the  Fourier  sequence.   Since  the  purpose  of  harmonic 
analysis  is  to  measure  the  contribution  of  all  frequencies  present,  another 
approach  must  be  considered. 

To  facilitate  a  more  through  approach  to  the  harmonic  analysis  reference 
is  made  to  an  analysis  by  U.S.  Carslaw  under  the  topic  of  practical  harmonic 
analysis  and  periodogram  analysis.   The  objection  sighted  above  was  overcome 
by  substituting  for  the  Fourier  series  a  trigonometric  series  with  a  limited 
number  of  terms.   This  is  done  in  the  following  manner.   Having  the  values 
of  the  time  series  for  one  period  given  at  the  points 
0,  a,  2a,  ...  ,  (m-l)a 
where  ma  =  211 
the  equidistant  points  on  the  time  axis  at  which  the  observations  are  taken 
are  denoted  by 
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V  *!»•••.  Vl 

and  the  corresponding  values  of  the  observations  by 

y0,  y±>  y2>  ■   •  •  >  vm_i 

The  time  series  is  then  represented  by  the  sequence 

f  (x)  =  a  +  a,  cos  x.  +  a„  cos  2x  +  .  .  .  +  a 
n       0    X      ±         &                1           n 

cos  nx 

+  b  sin  x  +  b  sin  2x  +  .  .  .  +  b 

sin  nx 

If  2n+l=m  the  fourier  coefficients  can  be  determined  sc 

i  that 

f(x  )  =  y     when  r  =  0,1,2,. ,.,2n 
r     r 

The  2n+l  equations  which  yield  this  determination  are 

a.  +  a,  +  .  . .        +a  +.  . .         +  a 
0    1             p                n 

=  yo 

a,  +  a,  cos  x,  + +  a  cos  pxn  +...  +  a  cos  nx  , 

0,1.    1        ,  p  .    1      ,  .  n  .    1} 
+  b.sin  x,  +...  +  b'sin  t>x,  +.  .  .  +  b  sin  nx, 
11        p    -  1         n      1 

=  yl 

a  +  a  cos  i  +...  +  a.  cos  px„  +...  +  a  cos  lix   , 

0.1.    2n     .  ,  p  .    2nt       n  .    2n) 
+  b.sin  x_  +...  +  Vsin  px„  +...  +  b  sin  nx. 
1     2n        p      2n        n      2n 

=  'y2n 

By  adding  the  above  equations 

2n 

(2n+l)a0=   J, 

r=0 

since  1+  cos  pa  +  cos  2pa  +...+  cos  2npa  =  0 

and   sin  pa  +  sin  2pa  +...+  sin  2npa  =  0 

when    (2n  +  l)  =  2H. 

Furthermore 

• 

1+  cos(pa)cos(ra)  +  cos(2pa)cos(2ra)  +. . . 

+  cos(2npa)cos(2nra)   =0          p  ^  r 

cos(pa)  sin(ra)  +  cos(2pa)  sin(2ra)  +...   p=l,2,. 
+  cos(2npa)sin(2nra)   =  0          r=l,2,. 

.  .n 
.  .n 
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and 


1  =  cos2(pa)  +  cos  (2pa)+...+  cos  (2npa)  =  5j(2n+l) 


Hence,  if  the  second  equation  is  multiplied  by  eoslpx^,  the  third  by 
cos(px  ),  etc.  and  added  the  result  will  be 

2n 
Js(2n+l)a  =  y  y  cos(pra) 

p   r^O  r 

In  a  similiar  manner  it  can  be  determined  that 

2n 
!-s(2n+l)b  =  I   y  sin(pra) 
P   r=l 

Hence,  a  trigonometric  series  is  formed  whose  sum  takes  the  required  values 
at  the  points 

0,  a,  2a,  ...,  2na  where  (2n+l)a  =  21! 
If  a  period  contains  an  even  number  of  observations  the  relationships  take 
the  following  forms .  The  interval  of  one  period  is  denoted  by 

0,  a,  2a,  ...,  (2n-l)a,  where  na  =  2n 
and  the  corresponding  values  of  the  observations  are 

yQ>  y1>  y2»---  y2n+i 

In  this  case  the  values  of  the  2n  constants  in  the  Fourier  series 

f  (x)  =  a  +  a,cos(x)  +  a„cos(2x)  +  ...  +  a   cos(n-l)x 
n       0    1         2  n-1 

+  a  cos(nx) 

n 

+  b,sin(x)  +  b„sin(2x)  +  ...  +  b   ,sin(n-l)x 
1         2  n-1 

So  that  this  series  yields  the  points  in  the  time  series  are  (9) 
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l     2n+l 

r=0 


,  2n-l 

=  „  I   yr  cos  Pra.  if  p  i  n 


p 


r=0 


2n-l 
"h   2n 


57  I  y„  cos  rn 


r=0 

2n-l 

bB  "  „  I     yr  sin  pra 

r=l 

where  a  =  n/n 

For  the  purpose  of  time  series  analysis  the  above  equations  while 
better  than  (l.lt.6)  are  still  inadequate.   The  expressions  above  depend  on 
the  criterion  that  the  number  of  observations  in  the  period  is  equal  to  the 
number  of  Fourier  coefficients  to  be  determined.   In  the  application  of 
Fourier  analysis  to  the  time  series  it  is  often  desirable  to  determine  the 
Fourier  coefficients  in  a  series  in  which  m  >  2n+l.   That  is,  the  number  of 
coefficients  in  the  series  is  less  than  the  number  of  observations  in  the 
basic  period.  There  are  two  justifications  for  the  above  condition 

(1)  The  computation  time  for  considerinc  all  m  frequencies  in  the 
Fourier  series  may  not  be  Justifiable. 

(2)  The  forecasting  technique  used  involves  matrix  manipulation. 
Since  the  calculation  time  for  matrix  manipulation  increases 
exponentially  with  the  order  of  the  matrix,  and  since  the  order 
of  the  matricies  increases  by  2  with  each  periodic  term  added  to 
the  forecasting  model,  these  terms  must  necessarily  be  limited. 
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For  the  purpose  of  determining  the  Fourier  coefficients  for  a  series 
in  which  in  >  2n+l  the  previous  analysis  vhich  depends  on  2n+l  =  m  is,  of 
course,  no  longer  applicable.   The  result  is  a  movement  from  the  deter- 
ministic evaluation  of  the  coefficients  to  an  approximation  of  the  coef- 
ficients based  on  the  theory  of  least  squares.  That  is,  the  coefficients 

V  V  •"'  V  V  ••"  bn 

are  determined  so  that  f  (x)  approximates  as  closely  as  possible 

yo>  yr  •■••  ym+i  at  V  V  •••'  Vi 

The  theory  of  least  squares  shows  that  the  closest  approximation  is 
obtained  by  making  the  function 

m-1 

I   (y  -f  (x  )f 
r=0  r  n  r 

a  minimum.  Where  the  above  sum  is  regarded  as  a  function  of 

V  V  •"'  V  \ bn 

The  conditions  which  make  this  sum  a  minimum  are  given  by  Carslav  (2) 

m-1 

I  (y   -f  (x  ))  =  o 


1 

£  (yr  "  W  C0S  pXr  =  ° 


r=0 


m-1 

£  (yr  "  fn(xr^  sin  pxr  =  ° 

r=0 


where  the  above  expressions  are  evaluated  for  p  =  1,  2,  ...  n.   The  con- 
ditions above  lead  to  the  following  values  for  the  coefficients 
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m-1 

maQ 

r=0  r 

m-1 

^p 

=  UT 

r=0  r 

m-1 

cos  prx 
P 

P 

r=l 

sin  prx 
P 

where  p  =  1, 

2,  ... 

n  and  m  is  odd. 

But 

if  m  is 

even,  the  coefficient  a  (where  p 
P 

=  5sn) 

LS 

m-1 

ma, 
55m 

=    lrr 

r=0 

cos  rll 

Although  these 

expressions  for  the  coefficients 

of  the  periodic  terms 

are 

adequate 

for  the  time  series  analysis,  one  more  improvement  can  be  made. 

Instead  of  limiting 

the  range  of  the  summations  from 

0  to  (m-l),  this  range 

can 

be  taken  as  the 

largest  integral  multiple 

of  the 

period  of  the  harmonic 

in 

the  time 

series. 

The  corresponding  equations  are 

given  by  Davis  (3,57) 

A(T)  ■ 

2  ?      2nt 

.  F  J^  yt  cos  — 

(1:U:7) 

B(T)  1 

N' 
2    r       ,   2Ilt 
■  F  JQ  vt  sin  — 

where  T  is  the  period 

of  the 

harmonic  considered 

N1  is  chosen  equal  to  the  largest  multiple 

of  T  in  the  time 

series 

• 
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1:5  The  Case  for  "Detrending"  the  Time  Series 

In  order  to  investigate  the  effect  of  the  secular  trend  of  the  mean 
on  the  harmonic  analysis  consider  the  case  of  a  linear  trend.  This 
development  is  given  by  Davis  (3,75). 

Suppose  the  time  series  in  the  interval 
-a  <  t  <  a 

has  the  trend 

y  =  yQ  +  mt 

Suppose  also  that  the  harmonic  analysis  reveals  that  the  time  series  has  a 
harmonic  term  of  the  form 

h(t)  =  A(T)  cos  2|i+  B(T)  sin^fi 

where  T  is  the  period  of  the  harmonic.  If  y,  above,  is  expanded  in  a 
Fourier  series  in  the  interval 

-  a  <  t  <  a 

the  result  is 


y  =  yo  +  ^<»inf-%sinfU|8in-|5*...,     (1:5:1) 

Now  if  in  h(t)  the  period  T  belongs  to  the  Fourier  sequence,  that  is,  if 
there  is  an  integer  n  such  that  n  =  2a/T,  then  the  corresponding  term  in 
(1:5:1)  must  have  been  included  in  the  coefficient  of  the  sine  term  B(T) 
obtained  by  the  Fourier  analysis.   Hence,  the  coefficient  of  sin(2nt/T) 
which  belongs  to  the  true  harmonic,  independent  of  the  trend,  must  be  B(T) 
diminished  by  that  part  due  to  the  trend. 
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Since  the  influence  of  the  trend  upon  the  harmonic  is  the  term 

,  . ,n  2ma  1    ,   ,n  mT 

;   n  n  "  l  ;   n 

the  true  harmonic  is  the  function 

h'(t)  =  A(T)  cos  S^+  B'(T)  sin  ^ 

where 

B'(T)  =  B(T)  +  (-l)n  f-                                                            (1:5. 

2) 

If  o2  is  the  variance  of  the  original  series,  then  the  variance  o2 

reduced  by  the  trend  and  the  harmonic  term  will  be 

°1  =  °2  "  °T  "  °H 

where  a2  is  the  variance  due  to  the  trend 

a2  is  the  variance  due  to  the  harmonic  term 

It  has  been  shown  in  (1 : U: U)  that 

o2  =  *s(A2(T)  +  B'2(T))  . 

For  the  trend 

n      2   3        J 

If  the  series  is  defined  over  the  interval 

0  <  t  <  2a 

instead  of  the  interval 

31. 


-  a  <_  t  _<  a 

then  the   only  modification  in  the   above   analysis   is   merely  that  B'(T-)    as 
given  in   (1:5:2)   is  replaced  by 

B'(T)   =  B(T)   +51 

In  the  case  in  which  2a/T  is  not  an  integer,  the  period  T  does  not 
belong  to  the  Fourier  sequence.  In  this  case  the  above  analysis  vill  yield 
only  an  approximation  to  the  reduced  variance  o^. 

The  above  analysis  for  a  linear  trend  can  be  easily  extended  to  other 
types  of  trends.  However,  trend  analysis  described  in  section  1:2  provides 
a  method  for  removing  the  secular  trend  of  the  mean  from  the  time  series . 
Hence,  if  the  data  is  in  the  "detrended"  form  then  harmonic  analysis  gives 
a  true  representation  of  the  periodic  component  of  the  time  series. 
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2.0  FORECASTING  THE  TIME  SERIES 

2.1  Moving  Averages 

Up  to  this  point  the  object  of  the  analysis  has  been  to  determine  the 
forecasting  model.   The  analysis  has  taken  the  following  forms 

(1)  Trend  analysis  using  regression  to  determine  the  trend 
of  the  mean  and  to  put  the  data  in  "detrended"  form. 

(2)  Autocorrelation  analysis  to  determine  the  basic  periocity 
of  the  time  series. 

(3)  Spectral  analysis  to  choose  the  periodic  terms  in  the  model 
The  next  step  in  the  analysis  is  forecasting  the  value  of  the  time 

series  for  the  future  period.  Considering  again  the  representation  of  the 
time  series 

y(t)  =  C(t)  +  e(t) 

where  y(t)  is  the  observed  value  of  the  time 
series  at  time  t. 

C(t)  is  the  process  which  the  time  series 
is  following. 
The  criterion  for  the  forecasting  technique  is  to  give  an  estimate  of 
the  process  by  effectively  damping  out  the  superimposed  noise.  Thus,  a 
technique  is  needed  which  will  seek  a  balance  between  the  ability  to  respond 
to  secular  changes  in  the  process  and  the  effect  of  error  in  the  forecast 
due  to  the  random  variation. 

A3  an  illustration  of  such  a  technique  the  moving  average  method  can 
be  considered.   In  this  case  the  process  is  considered  to  be,  at  least  locally, 
constant.   Hence,  the  process  can  be  described  locally  as 
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E(t)  =  a 
By  including  the  random  noise  the  time  series  can  be  represented  as 


The  technique  of  the  moving  average  gives  the  following  estimate  of  the 
process  (7,96) 

yt  +  yt+l  +  yt+2    '  '  '   yt+N-l 

Mt  =         5 

where  M  is  the  average  of  the  H  most  recent  periods  in  the  data.  If  the 
autocorrelation  analysis  (1.3)  shows  that  the  time  series  is  basically 
periodic  then  for  the  best  results  K  should  be  some  multiple  of  the  period 
in  order  to  negate  the  effect  of  periocity  on  the  value  M  . 

The  rate  of  response  of  the  moving  average  is  controlled  by  the  value 
of  N.  Since  each  of  the  N  most  recent  observations  is  given  the  weight 
1/H,  as  N  is  increased  the  response  of  the  model  to  the  most  recent  obser- 
vation is  decreased.  This  response  can  be  seen  more  clearly  by  writing  the 
recursive  relationship  for  the  moving  average. 


yt  ™  yt-N 


t    t-1      N 

Suppose  that  the  time  series  is  following  a  constant  process  with  super- 
imposed noise  about  a  mean  a' .   Then  suddenly  the  process  Jumps  to  a  new 
mean  a".   Brown  notes  that  it  would  take  N  observations  for  the  moving 
average  to  fully  adapt  to  this  change.  (1,99) 
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2.2  Exponential  Smoothing 

The  lag  in  the  rate  of  response  to  a  change  in  the  process  is  one  of 
the  most  critical  characteristics  of  the  forecasting  technique.  Moreover, 
even  though  the  rate  of  response  of  the  moving  average  can  be  altered  by- 
changing  the  value  of  H,  the  calculations  involving  this  change  must  be 
carried  out  over  the  vhole  range  of  the  observations.  Thus  a  more  palitable 
approach  to  the  forecasting  method  is  needed. 

The  recursive  relationship  for  the  moving  average 


yt  "  yt-N 


M  =  M    + 

t    t-1      H 

can  be  approximated  by  the  relation 


M 


t  "  1/N  *t  +  (1  "  1/N)  Mt-1 


where  M  is  used  to  denote  the  estimated  value  of  M  .  The  underying  as- 
sumption in  this  case  is  that  y    can  be  reasonably  estimated  by 
1/N  (M  . ) .  Brown  uses  S  for  smoothing  instead  of  M  for  moving  average 
and  obtains  the  relation  (1,107). 

St(y)  =  ayt  +  (l-o)St_1(y)  (2.2.1) 

where  s  (y)   is  termed  the  smoothed  statistic 
evaluated  at  time  t . 

a  is  an  undimensioned  ratio  similiar, 
but  not  equal  to  1/N. 
The  carrying  out  of  the  above  relationship  is  called  exponential  smoothing. 
By  rearranging  the  expression  for  exponential  smoothing  an  interesting 
observation  can  be  made.   Expressing  (2.2.1)  as 
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Bt(y)  =st_x(y)  ♦«(yt-8t_1(y)) 

the  current  value  of  the  smoothed  statistic  is  expressed  as  the  previous 
smoothed  value  plus  a  fraction  of  the  difference  between  the  value  of  the 
time  series  at  the  present  time  and  the  forecasted  value  at  the  previous 
period.   This  idea  of  updating  the  current  estimate  of  the  time  series  as 
a  function  of  the  error  of  the  previous  estimate  is  found  to  be  quite  con- 
sequential in  later  development. 

Up  to  this  point  no  sound  Justification  is  presented  for  the  use  of  the 
smoothed  statistic  as  a  representation  of  the  process.  This  Justification 
is  found  in  the  definition  of  expectation.  Hogg  and  Craig  define  expectation 
as  follows  (7).  Let  X  be  a  random  variable  having  a  P.D.F.  f(x),  and  let 
u(X)  be  a  function  of  X  such  that 

/"  u(x)f(x)dx 
exists,  if  X  is  a  continuous  type  of  random  variable,  or  such  that 
I   u(x)f(x) 

X 

exists,  if  X  is  a  discrete  type  of  random  variable.  The  integral,  or  the 
sum,  as  the  case  may  be,  is  called  the  mathematical  expectation  (or  expected 
value)  of  u(X)  and  is  denoted  by  e(u(X)].   That  is 

E(u(X))  =  I   u(x)f(x) 
x 

if  X  is  a  discrete  type  of  random  variable 

Brown  uses  expectation  to  obtain  the  following  proof  of  the  validity 

of  exponential  smoothing.  (1,101)  The  following  expansion  is  first 

performed 
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S(y)  =  uyt  +  (l-«)(«Vt_a  +  (1-a)  st_2(y>) 

=  oyt  +  o(l-o)  yt_x  +  (l-a)2(ayt_2  +  (l-a)  St_3(y)) 
=  ayt  +  a(l-a)  y    +  a(l-a)  yt_2  +  .  .  . 
♦  a(l-a)n  yt_n  +  .  .  .  +  (l-a)*  yQ 


t-1 


=  a  I      (l-a)k  yt_v  +  (l-a)\ 


k=0 


The  expected  value  of  the  above  expression  is  then 
E(s(y))  =  a  I   fob 


o     -* 


E(y)  01  X  8k  -  — S—  E(y)  =  E  (y) 
0      (1-8) 


where,  for  convenience,  B=(l-a).  Hence,  the  Justification  of  using  the 
smoothed  statistic  as  a  forecast  of  the  time  series  lies  in  the  fact  that 
the  expected  value  of  the  smoothed  statistic  is  the  time  series . 

In  the  comparison  of  the  technique  of  the  moving  average  with  exponential 
smoothing  it  is  important  to  consider  the  weights  given  to  the  observations. 
In  the  case  of  the  moving  average  the  M  most  recent  observations  are  given  a 
weight  of  1/N  while  all  other  observations  are  given  weight  zero.   In 
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exponential  smoothing  the  current  observation  is  given  a  weight  of  a  and 
the  weight  of  all  previous  observations  decreases  geometrica-ly  with  age. 

Moreover,  in  the  moving  average  technique  N  observations  must  be  car- 
ried "on  the  books"  at  all  times.  This  can  be  a  disadvantage  when  N  is 
large  and  when  a  large  number  of  time  series  are  being  considered.  Further- 
more, the  moving  average  assigns  no  weight  to  any  observation  beyond  the 
last  N  observations  even  though  the  contribution  of  these  older  terms  may 
be  significant.  In  contrast,  exponential  smoothing  carries  in  one  word  of 
data  all  the  history  of  the  time  series. 

On  the  topic  of  sensitivity,  it  is  a  simple  matter  to  change  the  value 
of  a  at  any  time  and  thus  alter  the  response  of  the  smoothed  statistic. 
Again  the  moving  average  technique  falls  short.  Although,  the  response  of  the 
moving  average  can  be  altered  by  changing  the  value  of  N  this  change  neces- 
sitates recomputation  throughout  the  whole  range  of  the  data.   On  the  basis 
of  these  comparisons  further  consideration  of  the  moving  average  technique 
is  ignored. 

2.3  General  Theorem  of  Exponential  Smoothing 

The  previous  analysis  assumes  that  the  process  can  be  adequately  rep- 
resented by  a  constant  model.   The  next  step  is  the  application  of  exponential 
smoothing  to  models  other  than  the  constant  model.   In  the  application,  of  expo- 
nential smoothing  to  the  constant  model  the  recursive  relationship  (2.2.1) 

a  =  St(y)  =  ayt  +  (l-a)  S^ly) 

where  the  process  is  £(t)  =  a 
gave  a  method  of  re-evaluating  the  estimate  of  the  coefficient  in  the  model 
with  each  observation.   If  this  technique  is  extended  to  more  complicated 
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models  then  a  method  must  be  determined  for  recursively  re-evaluating  the 
coefficients  in  those  models.   Brown  develops  the  extension  of  the  model 
to  an  n   degree  polynomial  (1,132).  In  this  case  the  process  is  repre- 
sented by 

£(t)  =  aQ  +  ajt  +  a2/2  (t2)  +  .  .  .  +  an/n:  (tn) 

The  Taylor  series  expansion  about  the  tth  observation  yields  an  esti- 
mate of  the  future  observations  as  follows 

where  y    is  the  k   derivative  evaluated 

at  time  t,  (in  this  case  t  is  taken  to 
be  the  current  value) 

y.   is  an  estimate  of  y 

_(k)  m&   I 

y*    dtk  '* 

t  is  the  forecast  interval 
Thus  the  Taylor  series  expansion  yields  the  following  estimate  for  the  next 
observation 

n     k  M 

n  ty,      n   k 

t+T   .'.   k!     ,1    k. 
k=0         k=0   K  • 

Hence,  the  forecast  is  in  terms  of  the  current  estimates  of  the  derivatives 
of  the  model.   These  derivatives  correspond  to  the  coefficients  that  are 
required.   The  Mediate  goal  then  is  to  estimate  these  derivatives  through 
the  technique  of  exponential  smoothing. 
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The  expression  for  the  smoothed  statistic  is   given  in   (2.2.1)   as 

St(y)  =  otyt  +   (1-a)   St_x(y) 

If  this  is  referred  to  as  single  smoothing  and  double  smoothing  is  defined 
as 

s]2)(y)  =  asW(y)  +  (1-a)  s|f)(y) 

Then  multiple  smoothing  of  order  k  can  he  defined  as 

StW(y)  =  sM)(y)  +  (1-a)  SW(y)  (2.3.1) 

The  fundamental  theorem  of  exponential  smoothing  given  in  Appendix  D 
states  that  if  the  observations  y  can  be  represented  by  the  model 

k   (kl 

t+T  k=o  k! 

then  the  general  smoothed  statistic  can  be  represented  as 

k 

%WW-  !   ("D^fer  I  jVi«*Jli  (2.3.2) 

The  significance  of  the  fundamental  theorem  to  the  goal  of  representing 
the  time  series  by  a  polynomial  is  as  follovs . 

(1)  The  general  smoothed  statistic  was  defined  in  (2.3.1). 

(2)  For  any  polynomial  of  degree  n,  by  (2.3.2)  n+1  smoothed  statistics 
can  be  written  in  terms  of  the  n+1  derivatives. 
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(3)  Using  the  n+1  simultaneous  equations  from  (2),  the  values  of 

the  derivatives  y^  '  can  be  solved  for  as  linear  combinations  of 
the  smoothed  data. 
Hence,  a  method  is  provided  for  recursively  estimating  the  values  of 
the  coefficients  in  the  polynomial  model.  Looking  at  the  computational  ef- 
fort, however,  for  each  observation  in  the  series  the  n+1  smoothed  statistics 
have  to  be  recalculated  and  the  n+1  simultaneous  equations  solved  for  the 
n+1  derivatives.   The  next  effort,  therefore,  is  to  simplify  these  operations. 

2.U  Matrix  Representation  of  Exponential  Smoothing 

From  Appendix  D  the  fundamental  theorem  of  exponential  smoothing  can 
be  expressed  in  matrix  form  as 

S  =  Ma  (2.4.1) 

where  S  is  the  nxl  vector  of  smoothed  statistics 
a  is  the  nxl  vector  of  coefficients 
M  is  a  nxp  matrix  with  elements  involving  infinite 
sums  of  powers  of  the  smoothing  constant  (where 
by  virtue  of  the  fundamental  theorem  n=p) 
With  the  expression  for  the  fundamental  theorem  expressed'  in  the  form 
(2.!*.1)  the  vector  of  coefficients  can  be  easily  solved  for  as  (1,137) 

a  =  S  M"1 

where  M   is  the  inverse  of  the  corresponding 
square  matrix 
This  type  of  recursive  relationship  is  quite  adaptive  to  computer  pro- 
gramming. 
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Brown  shows  that  Just  as  in  the  case  of  the  constant  model  the 
validity  of  the  recursive,  relationship  for  the  coefficients  in  the  poly- 
nomial model  can  be  proved  by  using  expectation  (1,138).  Suppose  that  the 
time  series  can  be  described  by  the  polynomial 

yt  =  a  +  bt  +  ct  +  .  .  .  +  gt"  +  e 

■  5(t)  +  et 

where  e  is  the  noise  sample 
£(t)  is  the  process 
Since  smoothing  is  a  linear  operation 

S(y)  =  S(0  +  8(et) 
But  by  definition 

e(s(e))  =  0 
and  it  was  shown  that  exponential  smoothing  yields  the  expected  value  of 
the  data  so  that 

E(s(y)J  =  S(5) 
Up  to  this  point  the  only  forecasting  model  that  has  been  developed 
is  the  polynomial.   The  reason  for  beginning  with  this  limitation  is,  of 
course,  that  the  fundamental  theorem  expresses  the  general  smoothed 
statistic  in  terms  of  the  derivatives 


h       dtk  '* 
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which  are  the  coefficients  of  the  polynomial  model.  The  polynomial  model 
is  by  no  means  adequate  in  representing  the  time  series .  To  be  an  effective 
representation  of  the  process  of  the  time  series  the  model  must  contain 
terms  which  express  both  the  trend  component  of  the  mean  and  the  seasonal 
component .  As  seen  in  the  section  on  representing  the  time  series ,  the 
seasonal  component  is  most  adequately  described  by  sine  and  cosine  terras. 
Hence,  a  transition  must  be  made  for  recursively  estimating  the  coefficients 
in  a  more  general  model. 

Suppose  that  the  process  can  be  represented  by 

5(t)  =  «1*'1(t)  +  a2f2(t)  +  .  .  .  +  anf„(t) 

=  I  a.f.(t)  (2. k. 2) 

fa    i  1 

i=l   x 

Where  the  functions  f. (t)  are  of  the  types 

(1)  Polynomials 

(2)  Trigonometric  functions 

(3)  Exponential  functions 
(h)     Emperical  functions 

In  some  cases  it  might  be  advantageous  to  use  fitting  functions  that 
are  emperical  such  as  the  number  of  building  contracts  let  6  years  ago. 
The  only  criterion  that  these  emperical  functions  are  required  to  meet  is 
that  their  value  be  known  both  at  the  time  the  forecast  is  made  and  at  the 
time  in  the  future  for  which  the  forecast  is  required.  However,  emperical 
functions  lead  to  computational  difficulties  far  beyond  those  of  the  other 
types.  Both  because  of  these  computational  difficulties  and  for  reasons 
of  interest  the  emperical  functions  are  avoided  here. 
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The  method  for  this  general  case  is  required  to  serve  tvo  objectives 

(1)  Provide  a  simple  iterative  proceedure  for  revising  the  estimates 
of  the  coefficients  in  the  forecasting  model 

(2)  Provide  a  means  for  discounting  the  weights  given  to  the  obser- 
vations according  to  a  time  scale. 

The  expression  for  the  forecast  is  given  by  Brown  (l,l6l). 

y(t+t)  -  a1(t)f1(t+x)  +  a2(t)f2(t+i)  +  .  .  . 

+  an(t)fn(t+t) 

n 
=  I     a  (t)f  (t+r)  (2.1*. 3) 

i=l 

The  residual  in  this  case  is  defined  as 

y(T-j)  -  y(T-j)  =  e(T-j)  {2.k.k) 

where  y(T-j)  stands  for  the  model  in  which  the 
coefficients  are  evaluated  with  all  the  data  through  time 
T  but  with  the  model  evaluated  J  periods  earlier 
Appendix  B  gives  the  expression  for  the  coefficients  that  minimize 
the  sum 

f    2   2 

till  '*  ^  (2-U"5) 

where  wt  is  the  weight  given  the  residual  at  time  t  as 

>'•?»?'  F"1  (2.U.6) 

where   W  is  a  T  x  T  matrix  in  which  W. .  is  the 

11 

square  root  of  the  weight  w  given  the 
residual  for  time  i,  and  all  off  diagonal 
elements  of  W  are  zero. 
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2  is  an  n  x  I  matrix  of  elements 

f.(t),  the  value  of  the  i   fitting 
function  at  time  t. 
F   is  the  n  x  n  symmetric  matrix 

T   9 
F  =  (*W)(*W)'  =  I     vf  f(t)  f'(t) 

t-1  Z 
If  the  data  is  discounted  as  in  the  case  of  exponential  smoothing 
then  the  weight  w  in  expressions  (2.U.5)  and  (2.U.6)  above  must  satisfy 
the  relationship  (1,163) 

v2   =8J 
T-j 

and  (2.1*. 5)  becomes 

T  n  .  2 

I   8J(y(T-j)  -  I   a(T)  f  (T-j))  (2.U.7) 

j=l  1=1 

The  F  matrix  becomes 

t-1 

^u^'  =  I     B  t,(X-3)   fv(T-j) 
w        1=0 

Hence  a  method  for  discounting  the  weights  given  to  the  observations  and 
recursively  updating  the  vector  of  coefficients  is  developed  for  the  gen- 
eral model. 

Since  the  major  part  of  the  calculations  involved  in  updating  the 
vector  of  coefficients  is  the  formation  of  the  F  matrix,  this  matrix  is 
chosen  for  further  consideration.  The  value  of  the  F  matrix  depends  on 
three  factors . 
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(1)  The  total  number  of  observations  in  the  time  series,  T. 

(2)  The  fitting  functions  contained  in  the  forecasting  model 

f.(t)   i  =  1,  2,  .  ■  ■  ,  n 

i 

(3)  The  weighting  function  v  vhich  in  this  case  is  given  by  the 
relationship 

w2   =  BJ 
T-J 

Hence,  the  F  matrix  does  not  in  any  way  depend  on  the  values  of  the 

observations  in  the  time  series.  Brown  (1,163)  uses  this  independence 

to  develop  a  recursive  relationship  for  the  F  matrix. 

F(t)  =  f(t)f'(t)  +  BF(t-l)  (2.U.8) 

Referring  to  expression  (2.4.6J  the  next  computational  effort  to  be 
considered  is  the  formation  of  the  data  vector  defined  as 


g(T) 


g2(T) 


g„(T) 


y  V  Mg 


,th 


Then  the  i   component  of  the  data  vector  can  be  written  as 


T-l 


gl(T)  =  I     SJ  y  (T-j)  f^T-j) 


(2.U.8) 


J=0 

Brown  (l,l6U)  develops  a  recursive  relation  for  this  vector  as 
g.(T)  =  y(T)  f.(T)  +  8  Ei(T-l) 

Hence,  from  (2.U.6),  after  n  observations  the  coefficients  can  be  estimated 
by 
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a'  =  a(T)  =  g(T)  F-1(T) 
Hence,  a  reasonable  forecast  for  the  general  model  is  (1,16k) 

y(T+t)  =  a'(T)f(T+T) 

n 
=  J  a. (T)f.(T+r) 

i=l  x    X 

Stopping  for  a  moment  to  evaluate  the  progress  of  the  forecasting 
development  the  following  is  noted. 

(1)  A  scheme  is  developed  for  applying  exponential  smoothing  to 
the  case  in  which  a  locally  constant  process  is  assumed. 

(2)  The  technique  of  exponential  smoothing  is  extended  from  the 
case  of  a  constant  model  to  a  general  polynomial  by  means  of 
the  definition  of  general  smoothing  and  the  general  theorem 
of  exponential  smoothing. 

(3)  Discounted  multiple  regression  is  introduced.  This  technique 
enables  the  further  extension  from  the  case  of  the  general 
polynomial  to  a  model  which  contains  both  polynomials  and  tran- 
sendendental  functions 

2.5  Computational  Considerations 

At  this  point  the  development  of  the  forecasting  model  is  complete. 
That  is,  the  time  series  under  consideration  can  be  adequately  represented 
by  a  model  composed  of  polynomials  and  transendentals .   Next,  consideration 
is  given  to  improving  the  computational  efficiency  of  the  forecasting  scheme. 

If  a  comparison  is  made  between  the  method  of  estimating  the  coef- 
ficients in  the  polynomial  model  and  the  method  using  discounted  multiple 
regression  an  important  difference  is  noted.   In  discounted  multiple  regres- 
sion the  coefficients  of  the  model  are  estimated  with  respect  to  a  fixed 


50 


time  origin.  On  the  other  hand,  in  the  case  of  the  polynomial  model  time 
is  measured  with  respect  to  the  most  recent  observation. 

Applying  the  concept  of  the  moving  time  origin  to  discounted  multiple 
regression  and  taking  t+1  =  T  observations  Brown  (1,168)  notes  that  the 
error  criterion  in  (2. It. 7)  becomes 

t   ,  2 

min  I     eJ(f'(-j)a(t)  -  y(t-j)) 
J=0 

This  is  the  same  as  the  error  criterion  given  in  (2.1*.T)  except  that  time 
is  counted  with  the  current  value  as  the  origin. 

If  the  same  change  in  the  time  origin  is  applied  to  the  data  vector 
in  expression  (2. k. 8)  the  result  is 

g.(t)  =  I     Bjf.(-J)y(t-J)  (2.U.8) 

1     J=0 

and  in  the  same  manner  as  the  development  for  the  fixed  time  origin  the 
expression  for  the  coefficients  that  minimize  the  error  is 

a'  =  y  W2*2F-1  (2. J*. 9) 

with  the  criterion  that  there  be  at  least  n  observations.  The  coefficients 
in  the  forecasting  equation  are  estimated  as  before  by 

a(t)  =  F_1(t)g(t)  (2.1*. 10) 

2.6  Recursive  Fitting  Functions 

With  certain  types  of  fitting  functions  the  value  of  the  vector  of 
fitting  functions  can  be  obtained  as  linear  combinations  of  the  value  of 
that  vector  at  the  previous  time  period.   In  the  cases  in  which  this  recursive 
relationship  holds  the  functions  are  said  to  have  a  fixed  transition  matrix. 
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That  is,  there  are  a  set  of  coefficients  L   which  do  not  depend  on  time 
such  that 

fx(t+l)  =  Lnfl(t)  ♦  L12f2(t)  ♦  .  .  .  ♦  Llnfn(t) 

f2(t+l)  =  L2lfl(t)  ♦  L22f2(t)  ♦  .  .  .  ♦  L2nfn(t) 


f  (t+1)  =  L  ,f,(t)  +  L  „f„(t)  +  .  .  .  +  L  f  (t)  (2.6.1) 

n        nl  1       n2  2  nn  n 

If  the  transition  matrix  is  represented  as  L  then  (2.6.1)  can  be 

represented  in  matrix  form  as 

f(t+l)  =  Lf(t)  (2.6.2) 

The  only  restriction  placed  on  this  transition  matrix  is  that  it  have  an 
inverse  L      .   The  fitting  functions  for  which  such  a  transition  matrix  exists 
are  the  polynomials,  exponentials  and  sinusoids.  Hence,  if  the  transition 
matrix  is  specified  along  with  the  vector  of  fitting  functions  at  time 
t=0  then  the  value  of  the  vector  of  fitting  functions  at  any  other  time  t 
can  be  determined  by  the  relation 

f(t)  =  I^ftO)  (2.6.3) 

Three  types  of  transition  matricies  used  in  combination  are  found  to  be 
quite  useful  for  the  time  series  considered.  Brown  (l,l65)  gives  the  trans- 
ition matrix  for  a  polynomial  as  an  n  x  n  matrix  with  ones  on  the  diagonal, 
ones  in  the  first  element  to  the  left  of  the  diagonal,  and  zeros  everywhere 
else.   For  example,  the  transition  matrix  and  initial  vector  of  fitting 
functions  for  the  cubic  model  are 
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10   0   0 

110  0 
0  110 
0   0   11 


f(0) 


Note:  In  the  case  of  the  polynomial  model 

the  coefficient  is  the  binomial  coefficient 

t! 


(t-k)I  k! 


If  the  fitting  functions  are  trigonometric  both  the  sine  and  cosine 
of  each  harmonic  must  be  included  (see  Appendix  A).  Thus  the  fitting 
functions  are 


f  (t)  =  sin  wt 


f  (t)  =  cos  cut 


Brovn  (l,l66)  gives  the  transition  matrix  and  initial  vector  of  fitting 
functions  as 

f(0)  = 

The  third  type  of  transition  matrix  is  for  the  case  in  which  growing 
amplitudes  and  shifting  phase  angles  are  included  in  the  periodic  terms 
(see  Appendix  A).   Suppose  the  example  above  is  expanded  to  include  the 
fitting  functions 


f  (t)  =  t  sin  uit 


f,  (t)  =  t  cos   tot 


Brown  (1,166)  gives  the  transition  matrix  and  initial  vector  of  fitting 
functions  in  this  case  as 
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cos   lo  sin  to  0                0 

-sin  to  cos   oj  0                0 

cos  oj  sin  w  cos  w   sin  uj 

-sin  u  cos  w  -sin  w   cos  to 


f(o)  = 


The  rule  for  generating  transition  matricies  for  forecasting  models 
whose  fitting  functions  are  linear  combinations  of  the  three  types  of  fit- 
ting functions  described  above  is: 

RULE  FOR  GENERATING  GRAND  TRANSITION  MATRIX 

Place  the  basic  submatricies  of  the  three  types  described  above  on  the 
main  diagonal  of  the  grand  transition  matrix  and  fill  in  the  required 
positions  with  zeros . 

For  example,  suppose  the  model  chosen  to  represent  the  time  series 
is  a  growing  sinusoidal  model  with  a  harmonic.  The  mathematical  expression 
for  the  model  is  then 


£(T+x)  =  (ax  +  agt)  +  (a3  +  a?t)  sin(2nt/12) 


+   (a,    +  afit)   cos(2nt/12)   +  a     sin(l»Ht/12) 
+  a„   cos ( knt/12) 


where  the  basic  period  is  12.     The  basic  submatricies,   then,   used  in 
building  the  grand  transition  matrix  must  represent 

(A)  A  polynomial  with  two  degrees  of  freedom. 

(B)  A  growing  sinusoid  with  frequency 

2nt/12 

(C)  A  harmonic  sinusoid  with  frequency 

Unt/12 
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The  schematic  representation  of  the  grand  transition  matrix  is  given  in 
Fig.  2.1. 


Fig.  2.1  Schematic  of  Grand  Transition  Matrix. 


where  the  O's  represent  the  required  zero  level  terms. 

Making  the  appropriate  substitutions  the  grand  transition  matrix 

is 
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Using  the  transition  matrix  and  the  nev  method  of  counting  time  the 
calculations  for  revising  the  estimates  of  the  coefficients  in  the  general 
forecasting  model  can  be  simplified.  The  expression  for  the  data  vector 
(2.1*. 8)  can  he  written  in  the  following  form  (1,169). 

t 
g(t)  =  y(t)f(0)  +  I   SJ  f(-j)  y(t-j)  (2.6.1.) 

J=l 

If  successive  values  of  the  vector  of  fitting  functions  are  generated  with 
the  transition  matrix  by  the  relation 

f(-j)  =  L_1f(-j+l)  (2-6-5) 

expression   (2.6.1*)   can  be  written  as 

t  . 

g(t)  =  y(t)f(0)  +     I     8J  IT1  f(-J+l)  y(t-j)  (2.6.6) 

By  changing  the  index  of  summation  by  the  relation  k=J-l,  the  recursive 
relation  for  the  data  vector  can  be  written  as 

g(t)  =  y(t)f(0)  +  S  L_1  g(t-l)  (2.6.7) 

In  the  above  expression  the  effect  of  the  new  method  of  counting  time  can 
be  seen.  In  expression  (2.6.7)  the  current  observation  is  weighted  by  the 
function  vector  f(0).  In  the  previous  method  of  counting  time,  with  a 
fixed  time  origin,  the  current  observation  is  weighted  by  the  function  vector 

f(t). 

Two  major  improvements  in  the  calculations  are  generated  by  changing 
the  time  origin  and  using  the  transition  matrix. 
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(1)  The  only  time  dependent  value  in  expression  (2.6.7)  is  the 
observation.  Hence,  the  components  of  the  data  vector  no 
longer  depend  on  absolute  time  and  can  be  tabulated  as  constants. 

(2)  A  simple  recursive  relationship  for  the  data  vector  is  developed. 
Turning  attention  to  the  F  matrix,  Brown  (1,170)  determines  the  re- 
cursive relationship  as 

F(t)  =  I   BJ  f(-j)f'(-j)  =  F(t-l)  +  B*  f(-t)f'(-t)  (2.6.8) 

J=0 

Even  with  these  simplifications,  however,  the  computations  have  not 
yet  significantly  decreased.   The  major  advantage  in  changing  from  a  fixed 
time  origin  to  a  moving  time  origin  is  found  in  the  following  property  of 
the  F  matrix.  In  the  cases  considered  the  fitting  functions  are  either 
trigonometric  functions  or  polynomials  and  3  is  less  than  one.  Under  these 
conditions  B*  tends  toward  zero  faster  than  f(-t)  can  grow  so  that  the  F 
matrix  reaches  a  steady  state  condition.   Hence,  F  inverse  can  be  determined 
in  its  final  form  for  any  set  of  fitting  functions  of  the  types  considered. 
The  expression  used  to  describe  this  convergence  is 
F(t)  =  F(t-l) 

Brown  (1,170)  notes  the  following  properties  of  this  convergence 
criterion.  If  a  fitting  function  is  used  which  takes  the  form  of  a  decreasing 
exponential  i.e. 

f(t)  =  e-at, 
the  F  matrix  will  reach  a  steady  state  only  if  past  data  is  discounted  at 
a  very  rapid  rate,  that  is  if 

B  <  e"2a  . 
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Moreover,  if  the  fastest  growing  function  in  the  model  is  t  ,  then  the 
number  of  periods  taken  for  this  convergence  is  given  approximately  by 

7  ♦  (5.l)n 

~  (LB)0"95 

Since  the  steady  state  conditions  are  assumed  to  have  "been  reached  the 
time  notation  is  dropped.   Moreover,  if  the  steady  state  condition  is  as- 
sumed the  conditions  for  F  to  have  an  inverse  will  necessarily  insue. 
Therefore  the  coefficients  in  the  model  can  be  estimated  by 

a(t)  =  F-1g(t) 
Therefore,  the  forecast  of  future  observations  is  given  by 

y(t+t)  =  a'(t)f(T) 

=  (F-1g(t))'f(T) 

=  g'(t)  F_1f(T) 
=  g'(t)c(T) 
where  T  is  the  forecast  period 

g'(t)  is  the  transpose  of  the  current  data  vector 
c(t)  is  a  column  vector  of  coefficients  that 
depend  only  on  the  values  of  the  fitting 
functions  at  time  t,  but  not  on  absolute 
time 

2.7  General  Exponential  Smoothing 

The  simplification  of  the  calculations  for  the  general  model  up  to 
this  point  depend  on  the  convergence  of  the  matrix  of  weighted  fitting 


58 


functions  F(t).     This   convergence  expressed  as 

F  =  F(») 

requires  two  conditions 

(1)  Successive  values  of  the  vector  of  fitting  functions  can  be 
generated  by  a  fixed  transition  matrix  L 

(2)  The  origin  of  time  is  taken  at  the  present 

Furthermore,  the  data  vector  can  be  defined  recursively  by  the  relation 

(2.6.1.) 

g(t)  =  y(t)f(0)  +  8L"1  g(t-l) 
Using  these  results  the  recursive  estimates  of  the  coefficients 

a(T)  =  (a^T),  a2(T),  .  .  .  ,  ^(T)) 

used  in  the   forecast  equation 

y(T+t)   =  a'(T)f(t) 

(2.7.1) 

=     I     a   (T)f   (t) 

i=l     x 

can  be  obtained. 

From  Appendix  D  the  minimum  discounted  squared  residual  sum  is  at- 
tained when 

F(T)a(T)  =  g(T) 

Furthermore,  when  F(T)  has  an  inverse  the  vector  of  coefficients  can  be 
expressed  as 


a(T)  =  F_1(T)  g(T) 
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and  if  the  convergence  criterion  is  met,  namely 


f,  (t)  <  S-t/2  for  all  i 


the  F  matrix  reaches  a  steady  state 

F  =  F(-)  =  I   BJ  f(-j)f'(-j) 
J=0 

Brown  (1,177)  substitutes  this  minimum  steady  state  solution  into  the 

recursive  relation  for  the  data  vector  (2.6. H).  The  result  is 

Fa(T)  =  y(T)  f(0)  +  SlT1  Fa(T-l) 
If  this  expression  is  premultiplied  by  F~ 

a(T)  =  y(T)F_1f(0)  +  8F"1L_1Fa(T-l)  (2.7.2) 

This  expression  can  be  analyzed  in  the  following  manner.  Defining  the 
time  independent  vector 
h  =  F-1f(0) 


and  the  time  independent  matrix 

H  =  BF"1L-1F 
expression  (2.7.1)  can  be  rewritten  as 

a(T)  =  hy(T)  +  Ha(T-l) 
Considering 

and  postmultiplying  the  definition  of  the  F  matrix  by  L'   L 
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J=0 

=  i(F-f(0)f'(0))L' 
p 

Hence 

H  =  3F_1L-1F  =  (l-F_1f(0)f' (0}f  (0))L' 
But  the  h  vector  is  defined  as 

h  =  F-1f(0) 
so 

H  =  L'  -  h(Lf(0))'  ■  V  -   hf'(l) 
Moreover,  (2.7.2)  can  "be  written  as 

a(T)  =  hy(T)  +  Ha(T  -  l) 

=  hy(T)  +  L'a(T  -  l)  -  hf ' (l)a(T  -  l) 
Since  f'(l)  a(T-l)  =  y(T-l)  is  the  forecast  of  what  the  observation  at 
time  T  will  "be,  as  of  the  data  received  through  tine  T-l.  The  above 
expression  may  be  written  as 

a(T)  =  L'a(T  -  l)  +  h(y(T)  -  y(T  -  1))  (2.7.3) 

and  hence  the  final  form  for  the  recursive  relation  for  the  estimates  of 
the  coefficients  in  the  general  forecasting  model  is  expressed  as  a 
function  of  the  value  of  the  vector  of  coefficients  at  the  previous  period 
and  the  error  of  the  forecast  made  at  the  previous  period. 

Exponential  smoothing  provides  an  estimate  of  the  coefficients  a(T) 
from  the  observations  y(0),  y(l),  .  .  .  ,  y(l').  If  there  is  no  noise  in 
the  data  and  if  the  fitting  functions  represent  the  process  then  if 
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exponential  smoothing  is  applied  over  a  long  enough  range  the  computed 
values  of  the  coefficients  vill  equal  the  true  values.  By  definition, 
hovever,  the  time  series  does  contain  noise.  The  representation  of  the 
time  series  is  given  as 

y(t)  =  6(t)  +  e(t) 

where  y(t)  is  the  observed  value  at  time  t 
£(t)  is  the  process 
e(t)  is  the  noise  in  the  t   observation. 

The  distribution  of  e(t)  has  the  properties 
EUt)  =  0  , 

E(eie  )  =  0  for  i  i  } 

=  o2  for  i  =  J, 

where  a2  is  the  variance  of  the  noise 
distribution. 
Exponential  smoothing  yields  a  forecast  whose  expected  value  is  shown  to  be 
the  process,  S(t).  Brown  (1,393)  also  uses  expectation  to  investigate  the 
estimates  of  the  coefficients. 

Suppose  the  process  can  be  exactly  represented  by  some  linear  combin- 
ation of  fitting  functions 

5  =  a 
where  the  vector  a  is  the  true  set  of  coefficients.  The  expected  value  of 
the  forecast  when  the  least  squares  criterion  is  used  is  shown  to  be 

B(y  )  =  af(t) 
Substituting  the  least  squares  estimates  of  the  coefficients  {2.k,9)   into 
the  expression  for  the  expected  value  of  the  coefficients  the  following  is 
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obtained. 

E  (a')  =  E(y  W^F-1} 

But  since  the 

expected  value  of  the  time  series  is 

=  5W2.?F"1 

=  a'^W^F"1 
=  a' 

the  process 

Hence,  the  expected  values  of  the  estimates  of  the 

coefficients  are  the  true 

coefficients. 

The  estimates  of  each  coefficient,  therefore, 

is  some  distribution 

whose  mean  is 

the  true  value  of  that  coefficient. 

Next,  the  variance  of 

the 

distribution  : 

Ls  considered.  The  correlation  coefficient  between  two  random 

variables  X. 

sind  X„  as  defined  previously  is 

E  (X1  -  u1)  (X2  -  Ug) 

p12  = 

01°2 

where  y  is  the  expected  value 

of  X 

y„  is  the  expected  value 

of  X2 

o  is  the  variance  of  X 

a     is  the  variance  of  X 

Hogg  and 

Craig  (T)  define  the  covariance  of  these  variables  as 

E^  -  »x)    (X2  -  u2)) 

Moreover 

,  the  variance  -  covariance  matrix  is 

defined  by  Hogg  and 
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Craig  (7)  as  follows.  Given  the  N  random  variables 

Xl'  V  •  •  •  '  h 

and  calling  the  variance  -  covariance  matrix  V,  the  elements  on  the  principal 
diagonal  of  V  are  respectively  the  variances  0. .  =  a?,  i  =  1,  2,  .  .  .  ,  N. 
The  elements  not  on  the  principal  diagonal  of  V  are  the  covariances 

o. ,  =  p. ,0.0, 

Exploiting  the  idea  of  the  variance  -  covariance  matrix,  the  corres- 
ponding matrix  for  the  variation  in  the  estimated  coefficients  a  caused  by 
the  noise  e  »  y  -  £  in  the  observed  data  is 

E(a  -a)(a  -  a)'  =  F_1  W'2   (x  -  5) ' (x  -  K)   W2,? '   F_1 

Since  the  noise  is   defined  to  have  no  serial  correlation  and  assuming  that 
all  the  noise  samples  have  the  same  variance  a2,  then   (1,393) 

E(x  -   %)<    (x  -   t)  =  la2 
and  the  covariance  matrix  for  the  coefficients  reduces  to 


(2.7.10 


F-^W2   (7W2)1    F-1o2   =  F"1  K  F"1  o2 
e  e 

where  K  =   (^W2)    (7W2)' 

The  K  matrix,  then,  can  be  represented  as 


K-  I     82J  f(-j)  f'(-j)  (2.7.5) 

J=0 

Moreover,  the  variance  -  covariance  matrix  for  the  coefficients  can  be 

expressed  in  terms  of  the  variance  of  the  noise  as 


6U 
Vo2  =  F"1  K  F_1o2  (2.T.6) 

£  E 

The  (i,j)  element  of  V   is  the  covariance  between  a.  and  a 
l  J  ^      J 

cov{a. ,  a, }  =  V.  o2 

vhere  a2  is  the  variance  of  the  (uncorrelated)  noise.  The  variance  of 
e 

the  i   coefficient  is 


var{a. }  =  cov{a. ,  a. }  =  V..o2 
1        11     il  e 


Hence,  the  elements  V   are  the  variances  of  the  coefficients  expressed  as 
a  multiple  of  the  variance  of  the  noise  o2. 
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3.0  COMPUTER  PROGRAMS 

The  analysis  presented  in  the  first  two  chapters,  trend  analysis, 
autocorrelation  analysis,  spectral  analysis  and  exponential  smoothing 
can  only  be  carried  out  practicably  through  the  use  of  automatic  computing 
equipment.  The  facility  available  for  this  study  is  the  I.B.M.  lM.0  system 
at  Kansas  State  University.   This  system  consists  of  the  following  equipment; 
an  I.B.M.  lUlO  computer  with  kOK  storage  capacity,  a  1U23  card  reader  and 
card  punch,  a  ll*22  printer,  an  I.B.M.  lltOl  computer  and  seven  T330  magnetic 
tape  drive  units . 

The  1410  system  is  internally  programmed  with  PR-155.   This  system  al- 
lows programming  in  either  Autocoder  or  Fortran.  In  this  study  Fortran  is 
used.  The  processor  occupies  10K  leaving  30K  for  the  compiled  program  and 
the  calculations.  This  limitation  of  storage  would  be  prohibitive  to  the 
application  of  the  preceeding  analysis  except  that  the  PR-155  system  allows 
phasing  of  the  program.  Phasing  consists  of  writing  the  program,  which 
itself  is  too  large  for  the  available  memory  capacity,  into  parts,  or  phases. 
These  phases  are  then  run  independently  and  anything  that  must  be  retained 
from  one  phase  for  following  phases  is  read  onto  a  "scratch"  file.  At  the 
completion  of  a  phase  the  processor  automatically  clears  core  and  loads  in 
the  next  phase. 

Phasing  the  program  does,  however,  have  its  associated  limitations. 
There  are,  basically,  two  limitations  which  must  be  considered.  The  first 
limitation  of  phasing  is  concerned  with  time.  Since  each  phase  is  compiled 
independently,  compiling  time  is  increased.  Moreover,  running  time  is  in- 
creased because  the  data  which  must  be  retained  between  phases  must  be  written 
on  and  read  off  tapes.  The  second,  and  most  crucial  limitation  is  concerned 
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with  the  nature  of  the  program.  Phasing  is  limited  to  programs  which  can 
he  devided  into  independent  parts  which  can  be  accomodated  in  the  memory 
capacity  of  the  available  facilities. 

The  programs  incorporated  in  this  analysis  provide  a  twofold  function. 
First,  these  programs  serve  the  function  of  a  medium  for  carrying  out  and 
evaluating  the  analysis  of  the  first  two  chapters.  Secondly,  these  programs 
serve  the  function  of  analyzing  the  practicability  and  economy  of  applying 
the  forecasting  techniques  to  a  digital  computing  system  with  a  Fortran  type 
of  processor. 

3.1  The  program  for  general  exponential  smoothing 

In  the  application  of  general  exponential  smoothing  to  computer  simu- 
lation a  general  program  is  required.   This  program  must  have  the  ability 
to  (l)  perform  general  exponential  smoothing  on  a  time  series,  (2)  change 
the  forecasting  model  and  (3)  vary  the  significant  parameters  in  the  fore- 
casting model.   In  the  last  requirement  these  parameters  are  taken  to  be, 
(a)  the  basic  period  of  the  model  and  (b)  the  value  of  the  smoothing  constant. 
Moreover,  it  is  required  that  the  program  for  the  application  of  general 
exponential  smoothing  perform  these  functions  in  a  reasonable  amount  of  time, 
with  the  ability  to  handle  a  wide  range  of  fitting  functions  and  time  series 
and  provide  output  in  the  desired  form. 

Considering  the  operations  involved  in  general  exponential  smoothing  the 
following  independent  phases  are  suggested.   The  first  phase,  phase  I,  is 
named  INCONT  and  executes  the  functions  of 

(1)  Reading  the  time  series  into  memory, 

(2)  Providing  the  description  of  each  forecasting  model  and 

(3)  Evaluating  estimates  for  the  initial  values  of  the  coefficients 
for  each  model  and  time  series  combination. 
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Upon  completion  of  the  first  phase  the  following  information  will  he 
available  to  the  insuing  phases 

(1)  Control  parameters 

(a)  The  basic  period  of  the  forecasting  model,  P 

(b)  The  number  of  fitting  functions,  H 

(c)  The  number  of  observations  in  the  time 
series,  ND 

(2)  The  time  series  X(l),  where  I  =  1,2,. ..,ND 

(3)  The  transition  matrix  TM(I,J) 

where  1=1, 2,..., H 

J=l,2 » 

(U)  The  initial  vector  of  fitting  functions, 
F(I),  where  1=1, 2,..., H 

(5)  The  initial  value  of  the  vector  of  coefficients  C(l), 
where  I  =1,2,...,N 

(6)  The  change  vector,  CHK(l),  1=1, 2,..., N. 

Evaluating  the  above  requirements  1,2  and  U  can  be  adequately  performed  by- 
normal  read  operations.  If  a  general  program  is  to  be  maintained,  however, 
requirement  3  and  5  must  be  considered  more  thoroughly.   Considering  re- 
quirement 3,  it  was  noted  in  Chapter  2  that  the  elements  of  the  transition 
matrix  which  are  included  to  describe  the  periodic  terms  in  the  model  take 
the  form 

+  cos  id  ,  +  sin  w 

where  oj  =  2II/P,  and  P  i3  the  basic  period  of  the  model.  Hence,  these  terms 
must  be  adaptable  to  a  change  in  the  basic  period  of  the  model.   Moreover, 
the  initial  values  of  the  coefficients  for  the  periodic  terms  in  the  fore- 
casting model  are  given  by  Brown  ( 1,19k)  as 


68 


P 

a  =  2/p  I  yk  sin  \ 
n      k=1  *     k 

for  the  coefficient  of  the  sine  terms  and 

P 

a  =  2/P  I   y  cos  u> 
n      k=1  k     x 

for  the  coefficients  of  the  cosine  terms.  Hence,  provision  for  calculating 
3  and  5  for  each  forecast  is  made. 

Considering  6,  the  change  vector  has  not  yet  been  discussed.  In  section 
2.6  the  recursive  relation  for  the  vector  of  fitting  functions  is  given  in 
(2.6.2)  as 

f(t+l)  =  Lf(t) 

where  L  is  the  transition  matrix  and  f(t)  is  the  vector  of  fitting  functions 
evaluated  at  time  t.  However,  because  of  the  moving  time  origin  the  F  matrix 
for  general  exponential  smoothing  is  defined  in  (2.6.1*)  as 

F(t)  =  I   Bjf(-J)f'(-J)  =  F(t-l)  +  BJf(-t)f'(-t) 

The  transition  between  the  recursive  relation  for  the  vector  of  fitting 
functions  defined  for  the  fixed  time  origin  and  the  recursive  relation  for 
the  same  vector  defined  for  a  moving  time  origin  results  in  expression  (2.6.5) 

f(-t)  =  L_1f(-t+l) 

Hence,  it  would  seem  that  the  requirements  of  the  F  matrix  necessitate  storing 
the  inverse  of  the  L  matrix.  The  objection  with  using  the  inverse  of  the 
transition  matrix  is  based  on  the  recursive  relation  for  the  vector  of  coe- 
fficients (2.7.3) 
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(T)  =  L'a(T-l)  +  h(y(T)  -  y(T-l))  , 


This  relationship  is  executed  in  the  last  phase  of  the  program.  Hence, 
the  original  transition  matrix  and  its  inverse  would  have  to  be  carried 
in  the  program.  This  dual  storage  is  an  obvious  burden,  if  not  a  lim- 
itation, on  the  program. 

Thus,  the  change  vector  is  introduced  to  avoid  the  calculation  and 
storage  of  the  inverse  of  the  transition  matrix.  Advantage  is  taken  of 
the  trigonometric  identities 

cos(-x)  =  cos(x) 
sin(-x)  =  -sin(x) 
and  the  relation 

t  for  n  even 


(-t)' 


-(t  )  for  n  odd 


The  fitting  functions  used  for  the  time  series  considered  are  either  simple 
povers  of  t  or  trigonometric  functions  or  multiples  of  these  functions.  In 
the  case  of  the  general  polynomial 

y  =  aQ  +  a^t)  +  a2(t2)  +  a3(t3)  +  .  .  .  +  an(tn) 


the  vector  of  fitting  functions  is 

(1) 
(-t) 
(-t2) 
(-t3) 

(-tn) 


1 
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As  an  illustration  of  a  more  complicated  model  the  example  in  section  2.6 

can  be  considered.  The  expression  for  this  model  is  given  as 

C(T+t)  =  (a  +  agt)  +  (a3  +  a?t)  sin(2JIt/12) 

+  (a,  +  a,t)  cos(2Ht/12)  +  a  sin(Unt/12) 

+  a„  cos(Unt/12) 

The  vector  of  fitting  functions  in  this  case  is 

1 

(-t) 

sin(u)  (-t)) 

cos(u)  (-t)) 

(-t)  sind^-t)) 

(-t)  cos(u  (-t)) 

sin(u_(-t) ) 

costoigt-'t) ) 

where  id  =  2Ilt/12  and  io„  =  Unt/12.  From  these  examples  it  can  be  seen  that 

f^-t)  =  +  fjU) 

vhere  f.(t)  is  the  value  of  the  i   fitting  function  evaluated  at  time  t. 

In  order  to  facilitate  this  relationship  in  the  computer  program  the  change 

vector  CHK(I),  where  I  =  1,2,...,N,  is  introduced.  This  vector  satisfies 

the  relationship 

=  0   for  f.(-t)  =  f.(t) 

CHK(l) 

=  1  for  f.(-t)  =-f.(t) 
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Thus,  the  transition  from  a  fixed  time  origin  to  a  moving  time  origin  is 
made  without  the  use  of  the  inverse  of  the  transition  matrix. 

The  functional  value  of  this  vector  can  he  seen  in  the  following 
example.  Suppose  the  transition  matrix  is  of  the  order  10  x  10.  If  a 
percision  of  ten  decimal  places  is  used  this  matrix,  or  its  inverse,  would 
require 

10  x  10  x  10  =  1000 
core  locations  for  storage.  In  contrast  the  storage  of  the  change  vector 
of  fixed  point  numbers  requires  only  10  core  locations. 

Finally,  all  the  required  information  from  the  first  phase  is  written 
onto  a  work  tape,  making  it  available  for  following  phases.  The  flow  dia- 
gram for  phase  I  is  shown  in  Fig.  3.1  and  the  program  is  shown  in  Appendix  E. 

The  second  phase,  named  RAY,  has  three  functions,  they  are 

(1)  Calculating  the  F  matrix, 

(2)  Calculating  the  K  matrix,  and 

(3)  Checking  for  convergence  of  the  F  matrix. 

The  expression  for  the  F  matrix  as  given  in  (2.6.8)  is 

t 
F(t)  =  I   BJf(-j)f'(-j)  =  F(t-l)  +  Btf(-t)f'(-t) 

As  shown  in  section  2.6  this  matrix  converges  under  specified  conditions. 
That  is,  F  =  F(»),  or  F(t)  ■  F(t-l).   The  convergence  criteria  used  in  this 
application  is  that  suggested  by  Brown  (1,   ) 


Fi1(t)  =  Fil(t-1}       6 
_iJ ,  JJ <  -if,-6 
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for  all  i  and  J.  Moreover,  the  K  matrix  is  given  in  (2.7.5)  as 

K  =  I   S2Jf(-j)f(-j) 

1=0 

Again  a  recursive  relation  for  this  matrix  can  be  formulated  as 

K(t)  =  K(t-l)  +  B2tf(-t)f'(-t)  . 

The  calculations  required  to  form  these  matricies  are  quite  extensive. 
Therefore,  it  is  necessary  to  reduce  the  calculation  time  for  these  oper- 
ations.  The  first  reduction  in  computing  time  can  be  attained  by  noting 
that  both  the  F  and  K  matricies  are  symmetric.  Moreover,  the  recursive 
relation  for  the  K  matrix  can  be  expressed  as 

K(t)  =  K(t-l)  +  (Btf(-t)f'(-t))st 
By  defining  the  term 

8tf(-t)f'(-t)  =  Z(t) 
the  two  recursive  relationships  can  be  written  as 

F(t)  =  f(t-l)  +  Z(t) 

K(t)  =  K(t-l)  +  8*Z(t) 

But  since  Z  is  also  a  symmetric  matrix  the  calculations  are  reduced  to 
computing  half  of  Z  for  each  iteration.  Moreover,  by  taking  advantage  of 
the  recursive  relationship 

the  calculations  are  further  reduced. 

The  check  for  convergence  requires  the  formation  of  the  quotients 
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Z.J(t)/F.J(t-l) 

for  i=l,2,...,n  j=l,2,...,n,  and  checking  to  see  if  these  terms  exceed 
10~  .  The  calculation  for  this  convergence  check  can  be  reduced  in  two 
ways 

(1)  By  noting  that  the  convergence  criteria  must  he  met  for  all 
the  elements  in  the  F  matrix,  the  program  can  be  written  to 
exit  from  this  check  routine  as  soon  as  it  finds  one  element 
which  has  not  converged. 

(2)  The  operations  involved  in  checking  for  convergence  are  quite 
extensive.   If  this  check  is  made  each  time  the  F  matrix  is 
updated  the  calculation  time  will  be  increased  considerably. 
Hence ,  the  program  can  be  written  to  make  this  check  at 
specified  intervals.  These  intervals  are  taken  to  be  each 

50   iteration. 
The  flow  diagram  for  phase  II  is  given  in  Fig.  3.2  and  3.3  and  the 
program  is  given  in  Appendix  E. 

The  third  phase,  named  MATINV,  performs  the  functions  of 

(1)  Taking  the  inverse  of  the  F  matrix, 

(2)  Forming  the  h  vector  and 

(3)  Determining  the  variance  of  the  coefficients. 

The  recursive  relation  for  the  coefficients  using  general  exponential 
smoothing  is  given  in  equation  (2.7.3)  as 

a(T)  =  L'a(T-l)  +  h[y(T)  -  y(T-l)) 
where 

a(T)  is  the  estimate  of  the  coefficients  at  time  T 
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L  is  the  transition  matrix, 
h  is  the  h  vector, 

(y(T)  -  y(T-l))  is  the  error  of  the  forecast  made  at  the 
previous  period. 
The  vector  h  is  defined  as 
h  =  F_1f(0) 

where  F   is  the  inverse  of  the  F  matrix  and  f(0)  is  the  vector  of  fitting 
functions  evaluated  at  time  t=0.  It  is  found  that  the  most  convenient  method 
for  performing  the  three  functions  required  by  this  phase  is  to  write  the 
program  for  determining  the  h  vector  as  a  subprogram  of  the  main  program 
for  finding  the  inverse  of  the  F  matrix  and  calculating  the  variance  of  the 
coefficients.  Again,  taking  consideration  of  the  memory  requirements  it  is 
noted  that  in  future  operations  the  F  matrix  is  not  needed.  That  is,  only 
the  inverse  of  the  F  matrix  will  be  required  after  this  phase.   Hence,  space 
can  be  conserved  in  core  by  replacing  the  F  matrix  by  its  inverse.  This  is 
equivalent  to  reading  the  inverse  of  the  F  matrix  over  the  F  matrix. 
In  section  2.7  the  variance-covariance  matrix  is  defined  as 

Vo2  =  F"1KF-1a2 

E  E 

where  o2  is  the  variance  of  the  noise  distribution.  By  using  this  matrix 

E 

the  variance  of  the   coefficients   expressed  as  multiples   of  the  variance   of 

the  noise  can  be   obtained  from 

var{a. }   =   cov{a. ,a. }   =  V     a2 
1  11  ii   e 

That   is,   the  elements  on  the  diagonal  of  the  covariance  matrix  provide  the 

required  information.     It  is   found,  however,  more   convenient  to  calculate 

the   entire  V  matrix  than  to  calculate   the   diagonal  elements   alone.      The   flow 
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diagram  for  phase  III  and  the  subprogram  for  calculating  the  h  vector 
are  shown  in  Fig.  3.1*  and  3.5,  and  the  programs  are  shown  in  Appendix  E. 

The  fourth  phase  is  the  phase  which  actually  makes  the  forecast, 
With  the  h  vector  calculated  all  the  information  required  for  this  phase 
is  available.  The  initial  estimates  of  the  coefficients  are  obtained  in 
phase  I.   With  these  starting  values  the  format  for  forecasting  the  time 
series  is 

(1)  Make  the  forecast  according  to  (2.7.1) 

y(T+x)  =  a'(T)f(-r) 

=  I   a  (T)f  (t) 

i=l 

where  y(T+i)  is  the  forecast  for  one  period 
in  the  future , 

a(T)    is  the  estimate  of  the  vector 

of  coefficients  made  at  the  present 
period  , 

f(-r)    is  the  vector  of  coefficients  eval- 
uated at  time  x  -   where  t  is  the 
forecast  period . 

(2)  Update  the  vector  of  coefficients  in  terms  of 

(a)  The  previous  vector  of  coefficients  and  the 
forecast  period, 

(b)  The  error  of  the  forecast  made  in  the  previous 
period. 

The  recursive  relation  for  the  vector  of  coefficients  is 
given  in  equation  (2.7.2) 
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a(T)  =  L'a(T-l)  +  hfy(T)  -  y(T-l)) 

where  a(T)  is  the  present  estimate  of  the  vector  of 
coefficients, 
a(T-l)  is  the  vector  of  coefficients  used  to  make 
the  forecast  in  the  previous  period, 
L'  is  the  transpose  of  the  transition  matrix, 
h  is  the  vector  of  constants  defined  in 
section  2.7, 
y(T)  is  the  present  value  of  the  time  series, 
y(T-l)  is  the  forecast  made  at  the  previous  period. 
(3)  Return  to  step  (l) 
The  two  auxiliary  functions  of  the  fourth  stage  are  to 

(1)  Calculate  the  sum  of  squares  of  errors  and 

(2)  Calculate  the  variance  of  the  forecasts. 

The  expression  for  the  variance  of  the  forecasts  is  given  in  section  h.k   as 

N 


A(yf 


yt' 

2  -  t=l 


>2 


F    N 

l*t 

t=l  * 

where  y  is  the  observation  at  time  t  and  y  is  the  forecast  made  for  time  t. 
The  flow  diagram  for  phase  four  is  given  in  Fig.  3.6  and  the  Fortran  program 
is  given  in  Appendix  E 

The  last  phase  of  the  program,  phase  V,  makes  a  plot  of 

(1)  The  time  series, 

(2)  The  forecasts  and 

(3)  The  absolute  error  of  the  forecast. 
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This  phase  is  named  "PLOTTER".  The  actual  plotting  is  done  in  a  sub- 
routine "PLOTS". 

It  is  found  that  a  serious  limitation  is  imposed  on  the  length  of  the 
time  series  that  can  he  accomodated  in  the  program  if  the  value  of  the 
forecast,  the  observations  and  the  errors  are  stored  in  memory  at  one  time. 
In  the  preceeding  phase  the  forecast  and  the  error  of  the  forecast  were 
calculated  for  each  period.  However,  if  these  sets  of  values  were  stored 
in  this  phase  the  length  of  the  time  series  that  could  he  used  would  he  too 
short  for  the  investigation.   Hence,  some  method  must  be  devised  to  make 
the  values  of  the  forecasts  and  the  forecast  errors  available  to  the  plot- 
ting phase  without  storing  them  in  the  previous  phase.  In  the  case  of  the 
forecast  this  problem  was  solved  in  the  following  manner.  Referring  to 
Fig.  3.6,  as  each  value  of  the  forecast  is  calculated  in  phase  IV  it  is 
written  onto  a  work  tape.  There  will  be  then,  at  the  completion  of  this 
phase,  N  forecasts  on  the  tape  for  each  model  used;  where  N  is  the  number 
of  observations  in  the  time  series.   In  the  plotting  phase  these  values 
can  be  read  in  and  transformed  into  the  subscript  variable  FCST(l),  where 
1=1, 2,..., N,  making  them  available  for  the  subroutine  "PLOTS".   In  the 
case  of  the  errors  this  difficulty  is  overcome  since  the  subroutine  "PLOTS" 
calculates  the  absolute  value  of  the  errors  independently  along  with  plotting 
the  three  values. 

This  subprogram  provides  a  vertical  plot  of  the  three  variables.  The 
horizontal  axis  on  which  the  three  variables  are  measured  is  limited  by  the 
ll*22  printer.  The  printer  is  capable  of  printing  133  characters  on  a  line. 
Hence,  the  scale  for  the  three  variables  must  be  transformed  to  an  integer 
scale  with  a  range  of  from  zero  to  133.  In  this  case  an  upper  limit  of  130 
was  actually  used.  The  vertical  scale  contains  one  line  for  each  period 


78 


in  the  time  series.  This  scale  does  not  have  an  upper  limit. 

The  plotting  function  is  carried  out  in  Fortran  IV  by  means  of  a  write 
statement.  For  each  period  in  the  time  series  the  statement 

WRITE(3,1)  (MP(L),(L=1,130)) 
is  executed.  This  statement  can  be  interpreted  as  writing,  by  means  of 
the  11*22  printer  (symbolic  unit  3),  by  the  accomodating  format  (7),  the 
subscripted  variable  MP(L)  which  ranges  from  1  to  130.  Now  if  for  each 
line  of  the  plot  L  takes  on  only  three  values  corresponding  to  the  ob- 
servation, the  forecast  and  the  error  for  that  period  on  the  integer  scale, 
then  the  required  plots  can  be  obtained.  One  of  the  characteristics  of 
Fortran  IV  is  that  a  variable  can  be  set  equal  to  an  alphebetic  or  special 
character.  Hence,  if  a  method  is  provided  to  set  each  of  the  three  values 
of  the  subscripted  variable  MP(L)  equal  to  the  corresponding  symbol  for  the 
plot,  then  the  requirements  for  a  plotting  routine  are  fulfilled.  The  flow 
diagrams  for  the  program  "PLOTTER"  and  the  subprogram  "PLOTS"  are  shown  in 
Fig.  3.7  and  3.8  and  the  programs  are  given  in  Appendix  E. 

The  second  phase  of  the  program  is  the  most  critical  with  respect  to 
time.  The  controlling  factors  in  this  phase  are  the  value  of  the  smoothing 
constant  and  the  size  of  the  matrix  being  handled.  The  calculation  time 
increases  factorialy  as  the  order  of  the  matrix. 

The  relationship  of  the  value  of  the  smoothing  constant  to  calculation 
time  can  be  demonstrated  a3  follows.  The  recursive  relation  for  the  elements 
of  the  F  matrix  is  given  in  equation  (2.6.5)  as 

F(t)  =  F(t-l)  +  Btf(-t)f'(-t) 
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Since  for  any  given  set  of  fitting  functions  the  rate  of  growth  of 
f  (-t)  is  fixed  the  rate  of  convergence  of  the  F  matrix  depends  on  the 
speed  at  which  B  goes  to  zero.   From  the  experience  of  this  investigation 
a  value  of  the  smoothing  constant  of  0.7  will  yield  convergence  for  the  F 
matrix  for  a  given  set  of  fitting  functions  five  times  as  fast  as  a  value 

of  0.9. 

In  this  section  the  flow  diagrams  for  the  separate  phases  of  the 
forecasting  program  are  shown.  It  is  noted  that  phasing  the  program  depends 
on  the  ability  of  the  system  to  retain  information  between  the  phases  by 
means  of  "work"  or  "scratch"  tapes.  Moreover,  it  is  noted  that  phasing 
also  depends  on  the  ability  to  write  the  program  in  independent  parts  or 
phases  that  can  be  run  separately.  In  this  investigation  one  more  limitation 
on  phasing  a  program  is  found.  Phasing  the  program  depends  on  the  ability 
to  write  the  program  into  phases  in  such  a  way  that  the  information  from 
any  phase  can  be  made  available  for  the  following  phase  or  phases  which 
require  it.  This  limitation  can  best  be  illustrated  in  the  context  of  the 
program  under  investigation.  Due  to  the  complexity  of  the  internal  read 
and  write  statements  between  phases  no  attempt  was  made  to  represent  them 
in  the  previous  flow  diagrams  of  this  section.  Figure  3.9  shows  the  data 
transfer  statements  for  the  "scratch"  files.   In  this  program  three  "scratch" 
files  are  required.  These  tapes  are  referred  to  in  the  listings  of  the 
forecasting  program  in  Appendix  E  as  symbolic  units  5,6  and  7.   It  is  found 
that  it  would  not  be  possible  to  U3e  any  less  than  three  tapes  since  the 
available  memory  capacity  of  the  system  would  not  allow  the  required  data 
transfer. 

The  read  and  write  statements  shown  in  Fig.  3.9  are  given  in  the  pro- 
gram using  free  style  formats.  Thus  each  of  the  read  and  write  statements 
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refer  to  a  physical  record  representing  one  model.  This  technique  is 
necessary  to  permit  a  multimodel  program.  This  type  of  multiphase  pro- 
gramming is  considered  to  be  a  significant  contribution  of  this  investi- 
gation. 

As  a  check  on  the  validity  and  accuracy  of  the  results  obtained  from 
the  program  for  general  exponential  smoothing  Table  3.10  is  constructed. 
This  table  gives  the  values  of  the  h  vector  and  the  variance  of  the  coef- 
ficients for  several  models.   For  each  model  these  values  are  obtained 
using  three  values  of  the  smoothing  constant  0.70,  0.90  and  0.95.  The 
results  are  shown  to  6  decimal  places  and  the  difference  betveen  those  ob- 
tained by  Brown  (l,l81*-193)  are  shown.  This  table  indicates  that  the 
results  obtained  are  in  essential  agreement  with  those  obtained  by  Brown. 
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3.2  The  Program  for  Calculating  the  Autocorrelation  Function  and 
"Detrending" 

In  section  1.3  the  autocovariance  is  defined  as 

T 

R  W  =  I     y,y,  k  /  (T-k+i) 

where  k  is   the  lag   for  which  the  autocovariance 
is   calculated 
T  is  the  range  of  the  series 
y     is  the  series  which  has  been  adjusted  to 
have  an  expected  value  of  zero 
The  normalized  form  of  the  autocovariance  is   defined  as 

p(k)  =  RxxW/Rxx(o)   • 

This  is  the  autocorrelation  coefficient.  The  set  of  values  for  the  auto- 
correlation coefficient  for  all  lags,  k  =  +1,  +2,  ...  ,  is  defined  as  the 
autocorrelation  function. 

This  expression  as  it  stands  does  not  readily  lend  itself  to  computer 
programming.  A  method  for  adapting  the  calculation  of  the  autocorrelation 
function  to  computer  programming  is  developed  "by  Raymond  W.  Southworth  (8). 
This  method  is  based  on  the  definition  of  the  autocorrelation  coefficient 
as  , 


N-k         N-k    N-k 

((N-k)  I   (y,)2  -  (  I  y.)2)  1=1        1=1 

i=l        i=l  x 


i=1  -  .-    i=1  -   i=1  —         N^k  a       N-k     „.l/2 


(3.2.1) 
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If  the  data  is  first  put  in  the  normalized  form  then  this  expression 
reduces  to  the  one  above. 

The  value  of  expression  (3.2.1)  is  found  in  its  adaption  to  computer 
programming.  First  the  following  sums  are  defined. 

N-k  N-k  2 

\  =  J/i+k  \  -  j/i 


N-k  H-l 

Fk  "  J/i  c* =  J^Vi* 


Then  it  is  noted  that  recursive  relationships  can  be  developed  for  the 
first  four  sums. 

Tk  "  Tk-1  "  yk 

Fk  =  Fk-1  "  yN-k+l 

Sk  "  Sk-1  "  yk 
2 
Gk  =  Gk-1  "  7N-k+l 

The  flow  diagram  for  the  computer  program  for  calculating  the  autocorrelation 

function  is  given  in  Fig.  3.11.  This  flow  diagram  illustrates  how  these 

recursive  relationships  reduce  the  calculations  extensively. 

This  program  also  removes  the  trend  of  the  mean  from  the  time  series. 

The  "detrending"  is  performed  in  a  subprogram  named  "TREND".   As  described 

in  section  1.2  this  operation  basically  consists  of  fitting  a  polynomial 

regression  model  chosen  to  represent  the  trend  of  the  mean,  to  the  data. 
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This  regression  curve  is  then  subtracted  from  the  data.  In  the  time  series 
considered  in  this  investigation  the  linear  regression  model  is  found  to 
be  adequate .  This  model  can  be  represented  as 
y(t)  =  a  +  b(t) 

where  y(t)  is  the  value  of  the  regression 

model  at  time  t  and  a  and  b  are  the 
coefficients  to  be  estimated  by  re- 
gression 
The  expression  for  the  constant  term  in  this  model  is  given  in  Appendix  B 

as 

N 
a  =  I  y(i)/N 
i=l 

where  N  is  the  number  of  observations  in  the  time  series.  The  expression 
for  the  coefficient  for  the  linear  term  is  given  as 


N  K 

I   U)y«  -  (I)  I  y. 

.   i=l i=l  x 


I   [if  -   2(1)  I   (l)  ♦  I   (3 
1=1  1=1     1=1 


Since  the  independent  variable  in  this  case  is  the  uniform  time  axis  this 
expression  can  be  reduced  to 


N  N     N 

I   (i)y  -  1/H  I   (i)  I   y 


N    „        N 
I    [if  -   1/N(  I    (i))' 
1=1  1=1' 


The  flow  diagrams  for  the  program  to  calculate  the  autocorrelation 
function  and  the  subprogram  to  "detrend"  the  time  series  are  given  in 
Fig.  3.11  and  3.12  respectively.   The  programs  are  given  in  Appendix  E. 
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3.3  The  Program  for  Calculating  the  Power  Spectrum 

The  basic  functions  of  this  program  are  calculating  the  coefficients 
of  the  harmonics  of  the  time  series  and  calculating  the  energy  term  for 
each  frequency  in  the  spectrum.  The  expressions  for  the  coefficients  are 
given  in  section  l.k   as 

N' 

A(T)  =  2/N'  I   y  cos(2nt/T) 

t=0  * 

N' 

B(T)  =  2/N'  I   y  sin(2nt/T) 

t=0  Z 

vhere  T  is  the  period  of  the  harmonic. 

N'  is  the  largest  multiple  of  T  in  the  series. 
In  section  l.k  the  energy  term  for  measuring  the  contribution  of  the 
frequencies  is  given  as 

E(T)  =  Sfill 
2o2 

where  R2(T)  =  A2(T)  +  B2(T), 

T  =  the  period  of  the  harmonic, 
o2  =  the  variance  of  the  data. 
For  convenience  the  value  of  R(T)  is  taken  as  a  measure  of  the  contribution 
of  the  frequencies.  This  is  compatable  with  the  spectral  analysis  found  in 
the  literature. 

The  input  to  this  program  is  the  detrended  data  from  the  program  in 
section  3.2.  The  flow  diagram  for  the  program  to  calculate  the  power 
spectrum  is  given  in  Fig.  3.13  and  the  program  is  given  in  Appendix  E. 
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A(I),  I  =  1,K  THE  INITIAL  V ALU  -  OF  THE  VECTOR  OF 
OF  THE  VECTOR  OF  FITTING  FUNCTIONS 
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PERFORM  THE  RECURSIVE  RELATION 
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K.,(T)  =  K.,(T-1)  +  1 


t„ 


Fig.  3.2  Flow  Diagram  for  Phase  "RAY'1 
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UPDATE  THE  VECTOE  OF 
FITTING  FUNCTIONS  FOE 
A  FIXED  TIME  ORIGIN 
f(T)  =  Lf(T-l) 


MAKE  THE  TRANSITION  TO  1 
MOVING  TIME  ORIGIN 

-f.(T)  if  CffiC(l)  =  I 

.(T)  = 

L     f.(T)  if  CHK(I)  »  0 


bW"1* 


GO  TO  A 


Fig.  3.2  continued 
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CHOOSE 
THE  NEXT 
ELEMENT 


COMPLETE  TH1 

f  matrix  by 
symmetry 


COMPLETE  THE 
K  MATRIX  m 
SYMMETRY 


Fig.  3.3  Convergence  Check  for  'KAY 
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READ  IN  THE  VALUE 
FOR  A  MODEL  FROM 
THE  WORK  TAPE 


CALCULATE  THE  IN- 
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MATRIX 
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Fig.  3.5  Flow  Diagram  for  Subroutine  "HViiU" 
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READ  IK  THE  MODEL 
FROM  THE  WORK  TAPE 
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UPDATE  THE  VECTOR 
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Fig.  3.6  Flov  Diagram  for  "FC3T" 
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CALCULATE  THE 
CAST  ERRORS 

FORE- 

(?) 

PLOT  THE  TIME 

SERIES  , 

THE  FORECASTS 

AND  THE 

FORECAST  ERRORS 

Fig.   3.7  Flow  diagram  for  "PLOTTER" 
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READ  IN  THE  CONTROL  STATEMENTS 

(1)  MAXIMUM  AND  MINIMUM  K/  IGE 
OF  THE  VARIABLES  Ti 
PLOTTED 

(2)  NUMBER  OF  FOISTS  TO  BE 
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(a)  FORECAST 
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( c )  ERROR 
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J_ 


CONVERT  THE  VALUES  OF  '-"-.". 
FORECAST,  THE  ERROR  AMD 
THE  OBSERVATION  FOR  PERIOD 
I  TO  THE  INTEGER  SCALE 


SET  THE  VARIABLES  MP( J), 
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HAVE 
ALL  THE  POINTS  BEEN, 

PLOTTED: 


Fig.  'j.ti     Flow  diagram  for  subroutine  "PLOTS' 
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READ 

H  =  THE  NO.  OF  OBSERVATIONS 

.T(MAX)  =  THE  MAXIMUM  RANGE  OF  THE  PERIODS 


READ  lu  THE  TIME  SERj.1; 
Y(I)   I  =  1,N 


'initialize 


T 


CALCULATE  Ti;E  VALUE  OF  N' 
THE  LARGEST  MULTIPLE  OF 
T  IN  THE  TIME  SERIES 


CALCULATE 

IV 
A(T)    =   2/I'J'      I     Y(I)    C0S(2irI/T) 

1=0 


CALCULATE 

rl' 
3(T)    =  2/1,"      I     Y(I)    sin(2   I/T) 
I   0 


CALCULATE 


R(T)    =    ((A(T))^      (BCD)"}'5 


INCREMEHT 

T  =  T  *   1 


NO 


Fig.    3.9     Program  to  Calculate   Rower  SJpectrum 
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PHASE  I   INCONT 

WRITE  ONTO  WORK  TAPE  7 

X(I) 

=  The  Time  Series 

ND 

=  The  Number  of  Observations 

NTOU 

=  The  Basic  Period 

N 

=  The  Number  of  Fitting  Functions 

F(I) 

=  The  Initial  Value  of  the  Vector  of  Fitting  Functions 

C(I) 

=  The  Initial  estimate  of  the  Vector  of  Coefficients 

MN 

=  The  Model  Number 

CHK(I) 

=  The  Change  Vector 

TM(I,J) 

=  The  Transition  Matrix 

PHASE  II   RAY 

READ  FROM  WORK  TAPE  7 

THE  PHYSICAL  RECORD  WRITTEN  IN  PHASE  I 

WRITE  ONTO  WORK  TAPE  6 

N 

=  The  Number  of  Fitting  Functions 

ND 

=  The  Number  of  Observations 

C(I) 

=  The  Initial  Values  of  the  Vector  of  Coefficients 

F1(I) 

=  The  Value  of  the  Vector  of  Fitting  Functions  Evaluated 

for  One  Period  From  the  Time  Origin 

X(I) 

=  The  Time  Series 

BETA 

=  The  Value  of  the  Smoothing  Constant 

MN 

=  The  Model  Number 

TM(I,J) 

=  The  Transition  Matrix 

WRITE  ONTO  WORK  TAPE  5 

H 

=  The  Number  of  Fitting  Functions 

F(I) 

=  The  Initial  Value  of  the  Vector  of  Fitting  Functions 

(BETA)   =  The  effective  Value  of  the  smoothing  Constant 

F(I,J) 

=  The  F  Matrix 

K(I,J) 

=  The  K  Matrix 

Fig.  3. 

10  Internal  Data  Transmission  Between  Phases 
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Fig.  3.10  CONTINUED 

PHASE  III   MATINV 

READ  FROM  WORK  TAPE  5 

THE  PHYSICAL  RECORD  WRITTEN  ON  PHASE  II 

WRITE  ONTO  WORK  TAPE  7 

F  =  The  Inverse  of  the  F  Matrix 

h  =  The  h  Vector 

PHASE  IV   FORECAST 

READ  FROM  WORK  TAPE  6 

THE  PHYSICAL  RECORD  WRITTEN  IN  PHASE  III 

READ  FROM  WORK  TAPE  J 

THE  PHYSICAL  RECORD  WRITTEN  IN  PHASE  III 

WRITE  ONTO  WORK  TAPE  5 

ND  =  The  Number  of  Observations 

X(l)  =  The  Time  Series 

PHASE  V   PLOTTER 

READ  FROM  WORK  TAPE  5 

THE  PHYSICAL  RECORD  WRITTEN  IN  PHASE  IV 
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A. 


READ 

H  =  THE  HO.  OF  OBSERVATIONS  IN  THE  TIKE  SERIES 

P  =  THE  UPPER  LIMIT  OF  THE  AUTOCORRELATION  FUNCTION 


READ  THE  TIME  SERIES 
Y(I)   I  =  1,N 


CALL  SUBROUTINE  TO 
"DETRSND"  THE  DATA 


GO  TO  ROUTINE  TO 
NORMALIZE  THE  DATA 


INI 

tial: 

ZE 

k  = 

0 

k 

=  O 

'"'k 

=  0 

I  = 

0 

CALCULATE  THE  SUMS 

1=1        K 

1   (Y(I))2  =  sv 

1=1 

SET 

Fk  "  Tk 

\  =  % 

1,  =  1 

CALCULATE  THE  MEAN 
OF  THE  DATA  Y 


CALCULATE  THE  SUM 
N 


>  (Yd)) 


2 


CALCULATE  THE 
VARIANCE 

Sv=(  I    (Y(I))2/N)3S 
"   1=1 


DIVIDE  THE  DATA 
3Y  THE  VARIANCE 

Y(I)  =  Y(i)/s 


SUBTRACT 

THE 

ME 

AN 

FROM 

THE 

DAT^ 

Y(I) 

-  t 

=  Y( 

I) 

Fit:.  3.11  Flow  Diagram  for  Program  to  Calculate  Autocorrelation 
Function 


0 

100 

APPLY  THE  RECURSIVE 
RELATIONS 

Tk  =  Vl-Y(k) 

Fk  ■  Fk-rY(M-k+1> 
sk  =  sk-l-  lY(k)) 

Gk  =  Gk-1~  (Y(;,'-k+]-)>2 

FOR  k=l,P 

CALCULATE 
N-k 
C  =   I     Y(I)Y(l+k) 
*    1=1 

FOR  k=l,P 

CALCULATE   FOR  k=l,P 
(S-k)C  -  T,  F 

((N-k)Gk-(Fk)^) 
x  ((N-k)Sk-(Tt.)2) 

Fig.  3.11  continued 
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CALCULATE  THE  CONSTANT  TERM 
IN  THE  REGRESSION  LINE 

a  =  I    (Y(D)AJ 
1=1 


CALCULATE  THE  LINEAR  COEFFICIENT 

S  $  K 

I  (I)Y(I)  -  1/K  I      (I)   V   Y(I) 

I  (i)2-  ( ?  (i))2 

i=i         i=i 


REMOVE  THE  TREND  FHOK 
THE  DATA 

Y(l)  =  Y(I)  -  (l)o 
FOR  1=1, N 


PRINT  THE  DATA  AND 
DETRENDED  DATA 


RETURN  TO  THE  MAIH 
PROGRAM 


Fi(».  3.12  Flow  Diaf-ran  for  Subroutine  "TKEJID" 
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U  APPLICATIONS 

The  application  of  the  analysis  in  the  first  three  chapters  is  carried 
out  on  two  time  series.  The  first  time  series  under  consideration  is  the 
international  airline  data  (Fig.  !».l)  and  the  second  time  series  is  the 
sheep  production  data  (Fig.  k.2).     The  application  takes  the  following 
form.   First,  "detrending"  is  performed  to  allow  the  data  to  be  viewed 
without  the  effects  of  the  secular  trend  of  the  mean.   Second,  autocor- 
relation analysis  is  performed  (l)  to  determine  if  there  is  significant 
evidence  that  the  time  series  is  generated  by  a  process  and  (2)  to  determine 
the  basic  period  of  the  time  series.  Third,  spectral  analysis  is  applied 
to  determine  the  contribution  of  the  harmonics  contained  within  the  basic 
period.   Fourth,  the  time  series  is  represented  by  a  model.   Finally, 
general  exponential  smoothing  is  performed. 

Within  the  context  of  the  application  of  general  exponential  smoothing 
and  the  preceeding  analysis  the  effectiveness  and  sensitivity  of  each 
technique  is  investigated.   In  the  case  of  trend  analysis  the  questions  to 
be  answered  are: 

(1)  Does  the  "detrended"  form  of  the  data  permit  more  effective 
analysis  of  the  time  series? 

(2)  Can  trend  analysis  be  used  to  effectively  determine  the  trend 
of  the  mean  of  the  time  series? 

(3)  Can  trend  analysis  be  used  to  provide  quantitative  estimates 
for  the  trends  of  the  periodic  component  of  the  time  series? 

As  an  aid  to  the  trend  analysis  a  Fortran  program  written  for  the  I.B.M. 
Il4l0  computer  is  provided  which  plots  the  "detrended"  data  to  make  the 
results  more  useable. 
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In  the  case  of  the  autocorrelation  analysis  the  application  to  the 
two  series  provides  answers  to  the  questions 

(1)  Does  the  autocorrelation  function  give  a  true  measure  of  the 
basic  period  of  the  process? 

(2)  Can  the  autocorrelation  function  be  used  as  a  test  of  signif- 
icance that  a  process  exists? 

In  the  case  of  the  international  airline  data  (Fig.  U.l)  it  is  auite 
evident  that  a  process  exists.   Moreover,  it  is  quite  evident  that  the  basic 
period  of  the  periodic  component  of  the  time  series  is  around  twelve  months. 
In  this  case  autocorrelation  analysis  can  be  verified  in  view  of  the  expected 
results.   The  sheep  data  (Fig.  k.Z),   however,  is  not  quite  as  obvious.   In 
this  case  the  trend  of  the  mean  seems  to  be  following  a  decreasing  function, 
but  there  is  no  immediate  indication  of  the  basic  period  of  the  periodic 
component  or,  as  a  matter  of  fact,  that  a  periodic  component  exists  at  all. 
In  contrast  to  the  international  airline  passenger  data  in  which  the  results 
of  the  autocorrelation  analysis  can  be  immediately  evaluated,  the  results  of 
the  analysis  in  the  case  of  the  sheep  data  can  be  evaluated  only  in  the  final 
stage  of  the  analysis,  the  forecasting  stage,  when  the  effects  of  the  choice 
of  the  basic  period  and  the  forecasting  model  can  be  tested. 

Using  spectral  analysis  it  is  determined  in  section  l.ll  that  the  contri- 
bution of  each  harmonic  within  the  basic  period  can  be  measured.   This  analysis, 
it  is  proposed,  will  lead  to  the  most  optimal  selection  of  the  periodic  terms 
to  include  in  the  forecasting  model.   However,  no  evidence  is  given  in  section 
l.ll  that  the  application  of  spectral  analysis  to  an  actual  time  series  will 
lead  to  the  obvious  distinction  between  the  contribution  of  the  harmonics. 
That  is,  it  is  not  yet  shown  that  the  application  of  spectral  analysis  yields 
results  which  justify  the  selection  of  a  limited  number  of  harmonics  to 
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adequately  describe  the  periodic  component  of  the  time  series. 

In  section  1.5  the  effect  of  a  linear  trend  in  the  mean  of  the  time 
series  on  spectral  analysis  is  considered.  This  analysis  provides  suf- 
ficient evidence  that  the  "detrended"  form  of  the  data  should  be  used  for 
the  spectral  analysis.  The  application  of  spectral  analysis  to  the  "de- 
trended  form  of  the  two  time  series  under  consideration  investigates  the 
effectiveness  of  this  method. 

The  final  point  to  be  investigated  in  spectral  analysis  is  its  relation 
to  autocorrelation  analysis.  In  section  1.3  it  is  shown  that  the  autocor- 
relation function  has  a  local  maximum  at  ok  ,  the  periods  of  the  harmonics 
which  describe  the  periodic  component  of  the  time  series.  Hence,  if  spectral 
analysis  is  carried  out  over  a  wide  enough  range  the  basic  period  of  the 
time  series  as  indicated  by  autocorrelation  analysis  should  agree  with  the 
first  harmonic  of  the  series  as  indicated  by  spectral  analysis 

In  summary,  then,  the  application  of  spectral  analysis  investigates 
the  following  questions : 

(1)  Can  spectral  analysis  be  effectively  used  to  determine  the 
periodic  terms  to  be  included  in  the  forecasting  model? 

(2)  Is  the  use  of  the  "detrended"  form  of  the  data  an  effective 
method  of  spectral  analysis? 

(3)  Do  the  results  of  spectral  analysis  verify  the  results  of  auto- 
correlation analysis? 

Continuing  to  the  final  3tage,  the  analysis  of  general  exponential 
smoothing  which  forms  the  contents  of  Chapter  2  presupposes  the  following: 
(1)  An  adequate  choice  of  the  fitting  function  which  make  up  the 
forecasting  model  can  be  determined. 
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(2)  The  basic  period  of  the  forecasting  model  can  be  adequately- 
determined. 

(3)  A  proper  choice  of  the  smoothing  constant  can  be  determined. 

In  the  application  of  general  exponential  smoothing,  trend  analysis,  auto- 
correlation analysis  and  spectral  analysis  can  be  effectively  used  to  help 
determine  (l)  and  (2).   It  is  veil  recognized,  however,  that  the  results 
obtained  by  these  measures  are  only  approximations.  Therefore,  the  effects 
of  errors  in  these  parameters  should  be  investigated.   This  investigation 
takes  the  form  of  sensitivity  analysis. 

Furthermore,  no  method  is  available  in  the  existing  literature  for 
determining  the  optimal  value  to  the  smoothing  constant.  As  a  matter  of 
fact,  no  evidence  is  available  to  substantiate  the  existence  of  a  singular 
optimal  value  of  this  constant.  This  research  does  not  attempt  to  find  this 
optimal.  It  is  recognized,  however,  that  a  local  value  of  the  smoothing 
constant  might  exist  for  each  particular  time  series  and  set  of  fitting 
functions.  Hence,  a  parametric  investigation  of  the  smoothing  constant  is 
carried  out  for  each  time  series. 

The  application  of  general  exponential  smoothing  then,  investigates 
the  following  questions: 

(1)  How  sensitive  is  general  exponential  smoothing  to  the  choice  of 
fitting  functions  used  in  the  forecasting  model? 

(2)  How  sensitive  is  general  exponential  smoothing  to  the  choice  of 
the  basic  period  in  the  forecasting  model? 

(3)  How  sensitive  is  general  exponential  smoothing  to  the  choice  of 
the  smoothing  constant? 
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U.l  The  Application  of  Trend  Analysis 

In  the  case  of  the  international  airline  passenger  data  (Fig.  k.l)  a 
linear  trend  seems  to  best  represent  the  mean  of  this  time  series.  The  flow 
diagram  for  the  computer  program  to  "detrend"  the  time  series  is  shown  in 
Fig.  3.12  and  the  program  itself  is  shown  in  Appendix  E.  The  results  of 
this  program  give  a  value  of  112.0  for  the  constant  term  in  the  regression 
model  and  a  value  of  2.6  for  the  coefficient  of  the  linear  term.   The  detrended 
data  is  given  in  Fig.  U.3.  Line  A-A  which  is  given  as  y  =  112.0  suggests 
that  the  "detrending"  has  successfully  removed  the  secular  trend  of  the  mean 
from  the  data. 

The  vertical  dashed  lines  in  Fig.  k.3  are  constructed  at  twelve  month 
intervals  in  order  to  indicate  the  basic  period  of  the  data.  The  most 
important  information  gleaned  from  the  trend  analysis  of  the  international 
airline  data,  however,  is  the  trend  of  the  periodic  component.  Line  B-B 
connects  the  peaks  of  the  periodic  component.  Since  the  secular  trend  of 
the  mean  has  been  removed  from  the  data  the  slope  of  the  line  B-B  indicates 
the  rate  of  growth  of  the  amplitude  of  the  periodic  component. 

Figure  U.lt  is  the  plot  of  the  "detrended"  sheep  data.  In  this  case  the 
linear  trend  is  again  assumed  and  the  results  yield  a  constant  term  of  2207.0 
and  a  value  of  -12.2  for  the  coefficient  of  the  linear  term  in  the  regres- 
sion model.   The  line  A'-A'  which  can  be  represented  as  y  =  2207.0  also 
indicates  that  the  linear  model  is  a  good  representation  of  the  trend  of  the 
mean.   The  "detrended"  form  of  the  sheep  data  does  not  obviously  yield  any 
further  information  at  this  point. 
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U.2  Application  of  Autocorrelation  Analysis 

The  flow  diagram  for  the  computer  program  for  autocorrelation  analysis 
is  given  in  Fig.  3.11  and  the  computer  program  is  given  in  Appendix  E.  The 
autocorrelation  function  for  the  international  airline  passenger  data  is 
given  in  Table  k,l   and  the  plot  of  this  function  is  given  in  Fig.  4. 5.  It 
can  be  seen  that  this  function  reaches  its  maximum  at  12  periods  and  multiples 
thereof.  Thus,  the  autocorrelation  analysis  bares  out  the  results  expected 
from  the  trend  analysis.  The  peak  of  the  autocorrelation  function  is  a 
maximum  at  12  months  and  decreases  at  multiples  of  this  period.   This  decrease, 
or  decay,  in  the  maximum  values  of  the  autocorrelation  function  is  due  to  the 
trend  in  the  amplitude  of  the  periodic  component  and  further  serves  to  verify 
the  effectiveness  of  autocorrelation  analysis. 

In  the  case  of  the  sheep  data  the  autocorrelation  function  is  given  in 
Table  h.2   and  the  plot  of  this  function  is  given  in  Fig  k.6.     This  function 
reaches  a  local  maximum  at  25  and  again  at  4o  months.  The  maximum  value  at 
Uo  months,  however,  is  much  more  predominent.   The  range  of  the  autocorrelation 
function  for  the  sheep  data  is  taken  as  50  months .  This  range  is  not  extended 
because  the  available  range  of  the  data  is  only  73  months. 

1».3  Application  of  Spectral  Analysis 

The  application  of  spectral  analysis  is  performed  on  the  "detrended" 
data.  Hence,  any  distortion  of  the  analysis  due  to  the  secular  trend  of  the 
mean  is  eliminated.   The  computer  program  for  the  application  of  this  analysis 
is  given  in  Appendix  E  and  the  flow  diagram  is  given  in  Fig.  3.13.   As  pro- 
posed earlier  the  purpose  of  spectral  analysis  is  to  determine  the  contrib- 
uting frequencies  in  the  time  series.   The  plot  of  the  power  spectrum  for  the 
international  airline  data  is  given  in  Fig.  1».7  and  the  power  spectrum  is 
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given  in  Table  U .  3 -   The  power  spectrum  for  this  series  has  a  clear  maximum 
at  both  6  and  12  months .  This  analysis  indicates  a  12  month  period  with  a 
harmonic.   Moreover,  since  one  of  the  purposes  of  spectral  analysis  is  to 
verify  the  results  of  autocorrelation  analysis  this  basic  period  of  12 
months  obtained  is  quite  significant.   The  clear  distinction  of  the  6  month 
harmonic  indicated  by  spectral  analysis  suggests  a  more  thorough  analysis 
of  the  periocity  of  the  time  series  using  this  method. 

In  the  case  of  the  sheep  data  another  interesting  result  is  obtained. 
Figure  h.8   and  Table  U . It  give  the  power  spectrum  and  the  plot  of  the  power 
spectrum  for  this  time  series.   The  power  spectrum  reaches  a  local  maximum 
at  18  months  and  a  maximum  at  36  months.   In  this  case,  however,  the  results 
do  not  agree  with  those  obtained  by  autocorrelation  analysis  which  indicates 
a  basic  period  of  40  months.  At  this  point  no  conclusion  can  be  drawn  about 
the  validity  of  the  two  techniques.  The  discussion  of  the  comparison  of 
the  two  methods  must  be  curtailed  until  the  effectiveness  of  the  forecasting 
models  can  be  considered. 

k.h     Application  of  Exponential  Smoothing 

The  flow  diagram  for  the  computer  program  used  in  exponential  smoothing 
is  given  in  Fig.  3.1  through  3.7  and  the  program  is  given  in  Appendix  E. 
Using  this  program  the  results  of  trend  analysis,  autocorrelation  analysis 
and  spectral  analysis  can  be  investigated.  Moreover,  the  effect  of  the 
value  of  the  smoothing  constant  on  the  forecast  can  be  investigated. 

The  application  of  general  exponential  smoothing  is  aimed  at  determining 
the  parameters  which  influence  th^  forecast  and  the  sensitivity  of  the  fore- 
cast to  those  parameters.   The  parameters  chosen  for  investigation  are 
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(1)  The  fitting  functions  used  to  describe  the  process 

(2)  The  basic  period  of  the  forecasting  model 

(3)  The  smoothing  constant 

In  the  comparison  of  forecasts  a  measure  of  the  effectiveness  of  each 
forecast  must  be  provided.   Brown  (1,393)  suggests  the  following  measure 
of  the  effectiveness  of  the  forecast 


I   (yt-yt)2/(T-n) 

t=l  l   t 

where  y  is  the  observation  at  time  t 
v  is  the  forecast  for  time  t 
T  is  the  total  number  of  observations  in  the 

time  series 
n  is  the  number  of  fitting  functions  in  the 
forecasting  model 
Upon  consideration  of  this  measure  of  effectivenss  proposed  by  Brown  the 
following  objection  is  incurred.  Although  the  above  measure  of  effectiveness 
is  useful  for  comparing  forecasts  of  the  same  time  series  it  cannot  adequately 
compare  forecasts  between  different  time  series.   The  reason  for  this  in- 
adequacy lies  in  the  fact  that  as  the  size  of  the  observations  in  the  time 
series  increases,  the  size  of  the  term  above  increases  even  though  the  errors 
may  not  be  proportionately  as  large.   Hence,  a  new  measure  of  the  effectiveness 
of  the  forecast  is  devised  as 

f  (yt-y//  I  yt  ^•u-l) 

t=i    t    *     t=i  t 

with  this  new  measure  of  effectiveness  the  forecasts  for  different  time 
series  can  be  compared  on  an  equal  basis. 
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The  application  of  general  exponential  smoothing  begins  with  the  con- 
struction of  the  forecasting  model  (see  section  1.1).  This  model  represents 
the  process  which  generates  the  time  series.   In  the  selection  of  the  fitting 
functions  used  to  represent  the  process  the  results  of  trend  analysis,  auto- 
correlation analysis  and  spectral  analysis  are  utilized.  Beginning  with  the 
international  airline  data  two  models  are  chosen.  The  first  model  repre- 
sents a  growing  sinusoid  with  a  basic  period  of  12  months  and  the  second 
model  represents  a  growing  sinusoid  with  a  basic  period  of  12  months  plus 
a  harmonic  at  6  months. 

The  mathematical  representation  of  the  first  model  is 


y(T+x)  =  (a  +a  t)  +  (a  +a  t)  sin(2nt/12) 


+  (aj+agt)  cos(2nt/12) 


where   the  terms : 


(a     +  at)     represent  the  linear  trend  of  the  mean 

a     sin(2Ilt/12)   +   a,    cos(2Ht/12)      represent  the 

12  month  periodic  component 
(at)    sin(2nt/12)   +    (agt)    cos(2nt/12)      represent 
the  growing  amplitude  of  the  periodic 
component, 
The  last  set  of  fitting  functions   also  give  the  model  the   ability  to  adapt 
to  shifting  phase  angles.      In  the  case  of  the  second  model  the  only  dif- 
ference in  the  fitting  functions  is  that  the  terms 

a,   sin(l»nt/12)   +  a„   cos(  liITt/12) 
7  o 

are  included  to  represent  the  harmonic  at  6  months. 


Ill 

The  initial  estimates  of  the  coefficients  for  the  terms  which  repre- 
sent the  linear  trend  of  the  mean  are  obtained  from  the  regression  analysis 
as  a^  =  112.0  and  a  =  2.6.  In  the  case  of  the  coefficients  a  and  a,-  of 
the  terms  which  represent  the  growth  in  the  amplitude  of  the  periodic  com- 
ponent trend  analysis  is  used.  In  section  U.l  these  coefficients  are  esti- 
mated through  the  use  of  Fig.  It. 3  as  1.9.  The  initial  value  of  the  coef- 
ficients of  the  periodic  terms  a  ,  a,  ,  a  and  a„  are  estimated  by  the  method 

proposed  by  Brown  (l,19H). 
P 
\   =  2/p  I  yv   sin  u(k) 
*      k=l  * 

for  the  coefficients  a  and  a7  and 

P 
\   •  2/P  I   y.  cos  u(k) 

X     k=l  * 

for  the  coefficients  a,  and  a„. 

The  purpose  of  this  investigation  is  to  determine  the  increase  in  ef- 
fectiveness of  the  forecast  due  to  the  inclusion  of  the  harmonic  term  as 
indicated  by  the  spectral  analysis,  and  the  effect  of  the  value  of  the 
smoothing  constant  on  the  forecasts.  The  results  of  this  application  are 
shown  in  Table  It. 5.  The  results  are  reported  according  to  the  forecasting 
model,  the  value  of  the  smoothing  constant,  8  and  the  basic  period,  P,  of 
the  forecasting  model.  The  measure  of  effectiveness  given  for  each  fore- 
cast is  the  variance  of  the  forecast  (U.it.l). 

In  this  application  the  effective  value  of  the  smoothing  constant  is 
reported.   That  is, 

Bn  =  8 

(effective) 

where  n  is  the  number  of  degrees  of  freedom  in  the  model.  Along  with  the 

results  In  Fig.  It.  5  the  other  results  obtained  for  each  forecast  are  given 
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in  Table  k.6.     This  table  gives  the  following  values  for  each  forecast 


EFFECTIVE  BETA 


F  INVERSE 


h  VECTOR 


VARIANCE  OF  THE  COEFFICIENTS 

The  advantage  of  tabulating  these  results  is  that  once  the  h  vector  is  cal- 
culated for  one  set  of  fitting  functions  and  a  particular  value  of  the 
smoothing  constant  then,  since  the  value  of  this  vector  is  independent  of 
the  time  series,  it  can  be  used  to  forecast  any  time  series  where  the  same 
fitting  functions  and  smoothing  constant  are  used. 

The  results  of  Table  U.5  clearly  indicate  that  the  value  of  the  smoothing 
constant  and  the  choice  of  the  fitting  functions  are  significant  parameters 
in  the  forecast.  In  the  case  of  this  data  a  value  for  the  smoothing  constant 
of  0.70  is  far  superior  to  a  value  of  0.90  for  either  of  the  two  models. 
However,  it  is  interesting  to  note  that  the  growing  sinusoidal  model  fares 
better  for  a  value  of  the  discount  factor  of  0.90  than  doe-  the  harmonic 
model.  When  the  discount  factor  is  reduced  to  0.70  the  situation  is  reversed. 
The  most  important  information  obtained  from  these  results  is  that  the  best 
forecast  is  obtained  using  the  harmonic  model.   This  varifies  the  results  of 
spectral  analysis. 
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The  effect  on  the  trace  of  the  forecast  for  a  change  in  the  smoothing 
constant  can  he  seen  in  Figs.  U.9  and  It. 10.  These  figures  clearly  indicate 
that  the  peaks  of  the  periodic  component  of  the  forecast  are  cut  down  as 
the  value  of  the  smoothing  constant  is  increased.  This  analysis  of  the 
smoothing  of  the  peaks  of  the  forecast  is  one  of  the  most  important  outcomes 
of  this  investigation. 

In  Fig.  I*. 9  the  ability  of  the  model  to  adjust  to  a  phase  change  can 
be  clearly  noted.  At  the  beginning  of  the  forecast  the  trace  is  clearly 
out  of  phase  vith  the  observations.  However,  at  the  end  of  the  time  series 
the  trace  is  exactly  in  phase. 

In  the  case  of  the  sheep  data  three  forecasting  models  are  used 

(1)  Linear 

(2)  Linear  plus  sinusoid 

(3)  Linear  plus  sinusoid  plus  harmonic 

Along  with  the  investigation  of  the  choice  of  fitting  functions ;  and  the 
value  of  the  smoothing  constant,  in  this  model  the  choice  of  the  basic  period 
is  also  investigated.  The  investigation  of  the  basic  period  is  noteworthy 
in  this  case  because  of  the  results  of  the  autocorrelation  analysis  and 
spectral  analysis  noted  in  section  U.3. 

The  mathematical  representation  of  the  models  used  in  forecasting  this 
time  series  are 

LINEAR 

C(T  +  t)  =  au  +  a2t 

LINEAR  WITH  SUPERIMPOSED  SINUSOID 

t(T  +  t)   =  a±  +  &gt  +  a     sin(2Itt/P)   +  a^   cos(2nt/P) 


Ill* 


LINEAR  MODEL  WITH  SUPERIMPOSED  SINUSOID  AND  HARMONIC 
t(T  +  t)  =  a  +  ftgt  +  a  sin(2nt/P)  +  a^  cos(2nt/P) 

+  ar  sin(Unt/P)  +  a,  cos(Unt/P) 
where  P  is  the  "basic  period  of  the  forecasting  model  and  t  is  the  forecast 
period.   In  contrast  to  the  models  used-  in  the  international  airline  pas- 
senger data  these  models  are  presented  in  a  more  general  manner  to  allow 
for  the  parametric  study  of  the  choice  of  the  basic  period.  Table  k.J   gives 
the  results  of  this  study. 

Just  as  in  the  case  of  the  forecasts  for  the  international  airline  pas- 
senger data  these  results  indicate  that  the  value  of  the  discount  factor  and 
the  choice  of  the  fitting  functions  are  significant  parameters  in  the  fore- 
cast.  In  this  case  the  effectiveness  of  all  the  models  used  increases  as 
the  discount  factor  decreases . 

The  most  interesting  study  in  this  time  series,  however,  is  the  study 
of  the  choice  of  the  basic  period.   The  significance  of  choosing  the  correct 
basic  period  can  be  seen  from  Table  U.7  which  indicates  that  if  the  basic 
period  is  chosen  incorrectly  at  12  months  then  the  linear  model  is  more  ef- 
fective. Hence,  an  incorrect  choice  of  the  basic  period  negates  the  effect 
of  the  periodic  terms  in  describing  the  time  series  and  even  makes  the  in- 
clusion of  these  terms,  at  the  cost  of  increasing  the  calculations,  detrimen- 
tal to  the  forecast. 

With  the  results  of  the  application  of  exponential  smoothing  available 
the  results  of  autocorrelation  analysis  and  spectral  analysis  on  the  shee-o 
data  can  be  reinvestigated.   Although  the  basic  period  of  36  given  by  spectral 
analysis  is  more  effective  than  the  result  of  Uo  given  by  autocorrelation 
analysis  the  harmonic  model  suggested  by  spectral  analysis  is  less  effective 
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than  the  basic  sinusoidal  model.  These  discrepancies,  however,  can  be 
partially  explained  by  the  range  of  the  available  data.   The  international 
airline  passenger  data  in  which  the  results  of  the  analysis  are  compatible 
contains  12  periods  of  data  while  the  sheep  data  contains  barely  2.   Hence, 
it  is  felt  that  the  results  of  this  study  indicate  one  more  criteria  for 
the  use  of  general  exponential  smoothing.   The  data  should  cover  enough 
basic  periods  to  allow  the  smoothing  technique  to  become  effective. 
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Table  It. 5  The  Results  of  the  Forecasts  of  the  International  Airline 
Passenger  Data  Using  General  Exponential  Smoothing 


MODEL 

GROWING  SINUSOID   P  =  12 

BETA 

VARIANCE  OF  FORECASTS 

0.70 

2.977 

0.90 

8.77* 

MODEL 

GROWING  SINUSOID  WITH  HARMONIC 

P  =  12 

0.70 

1.681 

0.90 

13.3»t3 
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Table  It. 7  The  Results  of  the  Forecasts  of  the  Sheep  Data  Using  General 
Exponential  Smoothing 


LINEAR 

MODEL 

BETA 

VARIANCE  OF  FORECASTS 

0.75 

9-025 

0.90 

11 .  Ul7 

LINEAR  1 

MODEL  WITH  SUPERIMPOSED  SINUSOID 

P  =  12 

0.75 

11.001* 

0.90 

13.212 

P  =  36 

0.70 

7.550 

0.75 

7.61*1* 

0.90 

9.352 

0.95 

12.295 

P  =  Ho 

0.70 

9-393 

LINEAR  MODEL  WITH  SUPERIMPOSED  SINUSOID  AND  HARMONIC 

P  =  36 

0.70 

10.196 

P  =  III) 

0.70 

1  . 
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DltfCUlMJIOK  OK  RciiULTij 

The   results   of  Luis   investigation   indicate  that  exponential  smoothing 
can  be   used  successfully   una  practicably   in   forecasting  time   series, 
nowever,    the  effectiveness   of  exponential  smoothing  must  be   qualified  by 
three   conditions 

(1)  The   correct   choice   of  the  parameters   of  the    forecasting  model, 

(2)  An   adequate   choice   of  the   value   of  the   smoothing  constant  and 

(3)  An   adequate   range  of  the   data. 

If  any  of  these  conditions  is  not  met,  exponential  smoothing  will  not  yield 
accurate  forecasts. 

Brown  presents  the  application  of  general  exponential  smoothing  as  a 
means  for  forecasting  a  time  series.   This  investigation  is  concerned  with 
developing  techniques  to  aid  in  the  selection  of  the  forecasting  model 
and  demonstrating  the  sensitivity  of  the  forecasts  obtained  by  exponential 
smoothing  to  the  parameters  of  the  forecast.   Trend  analysis  was  developed 
and  it  was  found  that  this  technique  is  quite  useful  in  selecting  the  terms 
in  the  model  which  describe  the  trend  in  the  mean  and  the  trends  in  the  per- 
iodic component.   In  the  time  series  considered  linear  trend  removal  was 
used.   However,  general  polynomial  trend  removal  was  developed  and  it  is  ex- 
pected tnat  this  extension  can  be  very  useful. 

Spectral  analysis  is  another  technique  used  for  determining  the  fore- 
casting model.   This  technique  is  compared  with  the  autocorrelation  analysis 
presented  oy  Brown  and  found  to  be  much  more  effective  for  selecting  the  per- 
iodic terms  in  the  forecasting  model.   The  combination  of  spectral  analysis 
with  "detrending"  to  remove  the  effect  of  the  trend  of  the  mean  is  an  orig- 
inal contribution  of  this  research  and  is  considered  one  of  the  most 
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important  developments.   In  the  case  of  the  sheep  data  it  was  shown  that  an 
incorrect  choice  of  the  basic  period  of  the  forecasting  model  negates  the 
effect  of  the  terms  which  describe  the  periodic  component  of  the  time  series 
and  even  makes  the  inclusion  of  these  terms  detrimental  to  the  forecast. 
When  the  correct  choice  of  the  basic  period  is  obtained  by  spectral  analy- 
sis, however,  these  terms  increase  the  accuracy  of  the  forecasts. 

The  comparison  of  the  forecasts  obtained  from  the  airline  data  and  the 
sheep  data  indicates  that  the  range  of  the  available  data  is  of  primary  im- 
portance in  the  application  of  general  exponential  smoothing.   Through  the 
use  of  spectral  analysis  and  autocorrelation  analysis  it  was  determined  that 
the  sheep  data  contained  barely  two  basic  periods  whereas  the  airline  data 
contained  twelve  periods.   Kence,  it  was  determined  that  the  more  accurate 
results  obtained  in  the  case  of  the  airline  data  were  due  to  the  requirement 
that  exponential  smoothing  must  be  carried  out  over  a  sufficient  range  of 
the  data  to  be  effective. 

The  effect  of  the  choice  of  the  smoothing  constant  on  the  forecasts 
can  be  seen  in  the  case  of  the  airline  data.   When  a  value  of  0.7  was  used 
the  forecasts  successfull  adjusted  to  the  trends  in  the  data.   On  the  other 
hand  a  value  of  0.9  for  this  constant  relulted  in  forecasts  which  degenerated 
in  accuracy  with  time.   This  investigation  did  not  propose  a  method  for 
choosing  the  value  of  the  smoothing  constant  but  merely  points  out  that  a 
local  optimum  for  a  particular  time  series  and  set  of  fitting  functions  does 
not  exist  and  that  a  study  of  the  choice  of  this  constant  merits  further 
investigation. 

This  investigation  did  not  make  any  assumption  about  the  normality  of 
the  distribution  of  the  forecast  errors.   If  the  normality  assumption  was 
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made  then  a  statistical  test  of  significance  such  as  the  F  test  could  be 
made  on  these  errors.  It  was  noted  that  the  normality  assumption,  if  made, 
would  have  to  be  preceeded  by  the  assumption  that  the  forecasting  model  is 
a  true  representation  of  the  time  series  and  this  assumption  was  considered 
invalid. 

The  results  of  the  programs  presented  demonstrant  that  the  analysis  of 
the  data,  trend  analysis,  autocorrelation  analysis  and  spectral  analysis 
can  be  successfully  and  economically  carried  out.   Moreover,  the  program 
presented  for  general  exponential  smoothing  shows  that  a  general  program 
which  can  accomodate  changes  in  the  forecasting  model,  the  basic  period 
of  the  model,  the  value  of  the  smoothing  constant  and  the  time  series  can 
be  developed  and  is  a  valuable  decision  aid  in  industrial  and  economic 
situations. 
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APPENDIX  A   FOURIER  ANALYSIS 

The  following  definition  of  Fourier  analysis  Is  given  by  Harold 
Thayer  Davis  (2,6l) 

"The  problem  of  Fourier  series  is  that  of  representing 
a  function,  either  continuous  or  with  a  finite  number  of 
finite  discontinuities  as  exhibited  by  a  set  of  discrete  data, 
by  means  of  a  series  of  fundamental  harmonics." 

The  above  statement  refers  to  a  harmonic  in  the  form: 

y  =  A  cos(2IIt/T)  +  B  sin(2Ilt/T) 
where  T  is  the  period  of  the  harmonic.   This  expression  can  be  written 
in  the  form 

y  =  A2  +  B2  cos(2nt/T  -  a ) 
where  a  is  the  lag  angle  given  by 

a  =  arc  tan  B/A 
An  example  illustrating  the  representation  of  a  harmonic  is  given  iri  Fig.  1A. 
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Fig.  1A  Representation  of  a  harmonic  term 
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y  =  6  cos(2nt/12)  +  8  sin(2nt/12) 

=  f'V  (6)2  +  (8)2  ]  cos  ((am/12)  -  a) 

where  T  is  the  period  of  the  harmonic , 
1/T  is  the  frequency, 

(A2  +  B2)**   is  the  amplitude, 

a   is  the  phase  angle  - 

The  purpose  of  defining  the  harmonic  term  lies  in  the  ability  to 
represent  a  function  by  a  series  of  these  harmonic  terms.  The  series 
composed  of  these  harmonics  is  called  the  Fourier  series.  The  value  of 
the  Fourier  series  lies  in  the  following  basic  theorem 

If  f(t)  is  a  single-valued  function  which  has  a 
derivative  throughout  the  interval  -a  <  t  <  a  except  for 
a  finite  number  of  points  at  which  it  has  finite  discontinu- 
ities, and  for  other  values  it  is  defined  by  the  equation 

f(t)  =  f(t  +  2a) 
then  f(t)  can  be  represented  by  the  series 

y  =  SjA  +  A  cos(nt/a)  +  A2  cos(2nt/a)  +  A  cos(3nt/a)  +  .  .  . 

+  B  sin(nt/a)  +  B  sin(2nt/a)  +  B  sin(3nt/a)  +  .  .  . 

This  series  is  the  Fourier  series.  The  coefficients  are  determined  from 
the  integrals 

A  =  1/a  /a  f(s)  cos(nns/a)ds  , 
n       -a 

B  =  1/a  fa   f(s)  sin(nns/a)ds 
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In  the  case  above  the  continuous  series  is  considered.  When 
analyzing  data  such  as  a  time  series  a  transformation  to  the  discrete 
case  is  in  order.  Moreover,  it  is  often  useful  to  convert  from  the  sym- 
metric form  shown  above  to  another,  more  suitable  form.  Since  in  a  time 
series  the  data  is  given  over  the  range 

0  <  t  <  2a 
the  transformation 

g(t)  =  f(t-a) 

is  made  and  the  integrals  are  rewritten  as 

2a 
A^  =  1/a  ^     g(t)   cos(nnt/a)dt  , 

2a 
B^  =  1/a  /     g(t)   sin(nHt/a)dt  . 

If  the  data  is   given  in  the  form 

tv   f2,   f3,    .    .    .    ,   fN     . 

Then  the  Fourier  coefficients  can  be  conviently  represented  in  the  form 
given  by  Davis  (2.610. 

B 

An  =  2/N  I  ft  cos(2nt/N) 


B 
B  =  2/N  I     f  sin(2nt/W) 
t=l  t 
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APPENDIX  B  REGRESSIOH  ANALYSIS 

In  the  case  of  linear  regression  a  straight  line  is  fitted  to 

the 

data 

in  such  a  way  as  to  minimize  the  sum  of  the  squares 

of  errors  (6) . 

In 

this 

case  the  error  is  defined  as  the  difference  between 

the  actual  data 

.  point 

and  the  value  obtained  by  the  straight  line  at  that 

point. 

The  equation  for  the  straight  line  is  utilized 

in  the  form 

y'  =  a  +  b(x  -  x) 

B- 

-1 

where  b  is  the  slope  of  the  line  and  a  is  the  Y  int 

ercept  on  the  line 

x  =  x.  The  problem  is  to  determine  the  parameters 

a  and  b  so  that 

the 

sum 

of  the  squares  of  the  errors  of  estimation  will  be 

a  minimum.  Let 

the 

co- 

ordinates  of  the  i   point  be  denoted  by  (xi,yi). 

Then  the  term  to  be 

minimized  is 

I  br4  -  r[f 

i=l 

where  y!  is  determined  by  B-l. 
1 

If  this  function  is  denoted  by  G(a,b)  and  written  as 

G(a,b)  =  I     [y,  -  a  -  b(x  -x)j 

1=1  X 

then  the  conditions  for  the  above  expression  to  be 

a  minimum  are  that 

its 

partial  derivatives  vanish.  Hence  a  and  b  must  satisfy  the  equati 

ons 

|£  «  l   2  (y  -  a  -  b(x-x))  (  -l)  =  0 

|2.  -  I   2  (y  -  a  -  b(x-x))  [  -  (x-x))  = 

0 
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When  the  summations  are  performed  term  by  term  and  the  sums  that  involve 
y  are  transposed,  these  equations  assume  the  form 

an  +  b£  (x-x)  =  ][y 

&l   (x-x)  +  b£  (x-x)   =  £(x-x)y 
Since  £(x-x)  =  0,  the  solution  of  these  equations  is  given  by 

I(x-x)y 

a  =  y  and  b  =  ;r 

I(x-x)2 

For  computational  purposes,  it  is  convenient  to  change  the  form  of  the 
expression  for  b  in  the  following  manner 


£xy  -  x£y 
b  = 


lx2   -  2x£x  +  Ix 


£xy  -  nxy 


"  v  2    -2 
lx     -  nx 

The  concept  of  linear  regression  can  be  easily  extended  to  polynomial 
regression.  (6)  Let  the  degree  of  the  polynomial  be  k  and  let  the  equation 
of  the  polynomial  be  written  in  the  form 

*'    =  C0  +  C1X  +  C2X2  +  •  •  •  +  c/ 

Aa  in  the  case  of  linear  regression  the  unknown  coefficients  are  estimated 
by  the  method  of  least  squares.   This  is  equivalent  -to  minimizing. the  sum 
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I   (y.  -  y:t 

. L.   1    1 


Since  it  is  more  convenient  to  work  with  variables  measured  froir.  their 
sample  means  than  with  the  variables  themselves,  the  following  definitions 
are  made 

y  =  Y  -  Y 

x.  =  X.  -  X. 

J     J     J 

If  y'  is  defined  by  y'  =  Y'  -  Y,  then 

Y  -  Y'  =  y  +  Y  -  (y'  +  Y)  =  y  -  y' 
If  now  the  capital  X's  and  Y's  are  expressed  in  terns  of  the  small  x's  and 
y's  the  polynomial  regression  model  can  be  written  in  the  form 

y'  =  aQ  +  a1x  +  a2x  +  .  .  .  +  a,  X  £-2 

Since  minimizing  £(Y  -  Y' )      is   equivalent  to  minimizing   £(y  -  y')      it   is 
just  as  well  to  determine   the   a's   to  minimize  the   latter  sum  which  because 
-of  13-2  may  be  written  as 


G(aQ)    a^    .    .    .    aj    =   J[  (y  -   aQ  -   <yc  -    .    .    .    -   a^) 


2 


If  this  function  is  to  have  a  minimum  value,  it  is  necessary  that  its  partial 
derivatives  vanish  there.  Hence,  the  a's  must  satisfy  the  equations 


3G  _  _3G  _        _3G 
3aQ   3&1  "  •  •  ■  "  3^  " 


This  differentiation  produces  the  equations 


157 


l2{y   -  aQ  -  axx  .  .  .  Vk](-l)  =  0 
My  -   aQ  -  alX-  .  .  .  -  a^)  (-x)  =  0 

[2(y-a0-  V-  .  .  .  -  Vk)[-J)    =0 

If  these  equations  are  multiplied  by  \,   the  summations  performed  term  by- 
term,  and  the  first  sum  transfered  to  the  right  side,  these  equations 
will  assume  the  form 

V  +  aJx  +  •  •  •  +  \Ix   ■  b 

a0Ix  +  allx     +   •    •    •   +  \lx         ■  [xy 


V  k  r  k+1  ^  r  2k        V  k 

a0Zx     +  ax2.x         +   •    •    •   +  a^x       =  £x  y 


Since 


and 


5>J  =  £(XJ  -  X)  =  0 


^  -  Id  - 1)  ■  o , 

all  terms  in  the  first  equation  except  the  first  term  vanish.  This  implies 
that  a  =  0,  and  thus  the  number  of  equations  to  be  solved  has  been  reduced 
by  one.  The  problem  is  now  reduced  to  solving  the  equations 
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ajx2  ♦  a2[x3  ♦   .    .    .  +  ajxk+1  =  J 


xy 


V       a2^-x      •  •  •     Vx     =  lx  y 


r  k+1     r  k+2  .  r  2k   r  k 

\l   x    +  a2ix    +  •  •  •  +  ajx   =jxy 

Having  developed  the  expressions  for  the  coefficients  in  the  regression 
model  for  the  general  polynomial,  the  next  step  is  to  consider  a  more 
general  regression  model.  Suppose  the  model  contains  terms  of  the  following 
types 

(1)  polynomials 

(2)  trigonometric  functions 

(3)  exponential  functions 
(U)  emperical  functions 

In  this  case  the  vector  of  coefficients  can  be  represented  as 


and  the  fitting  functions  are  represented  in  vector  notation  as 


f^t) 


f(t) 


f2(t) 


tjt) 


159 


where  f.(t)  is  the  value  of  the  i   fitting  function  evaluated  at  time  t. 
Hence,  the  data  is  represented  by  the  model 


X.  »  I     a.f.(t)  =  a'f(t) 
t   , fc,   1  i 

1=1 


The  criterion  for  the  selection  of  the  coefficients  is  taken  as 

T 

V   2  2 

mln  L     wt  et 

t=l  *  ' 

where  e(t)  is  the  residual  defined  as 

e(t)  =  x(t)  -  a'f(t) 

and  w  is  the  weight  given  the  residual  at  time  t. 

The  above  expressions  can  easily  be  put  into  matrix  form.   A  matrix 

—r  i      \  th 

? is  defined  as  an  m  x  T  matrix  of  elements  f. (t),  the  value  of  the  i 
fitting  function  at  time  t.  A  row  vector  x  is  defined  to  be  the  sequence 
(x.  ,  x  ,  .  .  .  ,  x,^)  of  values  given  by  the  model  a'   .  Moreover,  e  is 
defined  as  the  sequence  (e.,  e  ,  ...  ,  e  )  of  residuals,  where 

et  =  xt"  \   =  xt"  »•<*>«*) 

In  order  to  find  the  expression  for  the  coefficient  vector  that 
satisfies  the  regression  criteria  matrix  notation  is  continued.   Let  V  be 
a  TxT  matrix  in  which  W. .  is  the  square  root  of  the  veight  given  the 
residual  at  time  i.   All  off  diagonal  elements  of  W  are  zero-   Noting  that 
the  expression  of  the  model  is 

a' y  m   x  -  e 

the  residuals  can  be  expressed  as 
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e  =  x  -  a1  2". 

New  for  a  particular  choice  of  a  the  sum  of  squares  is 
Sa  =  ee'  =  (xW  -  a',?W)(xW  -  a'^W)' 

=  x2  WW  -  2a'xS'WW'  +  (a')2  (S«)(3V) 
The  particular  set  of  coefficients  that  minimize  this  sum  is  found  from 


oa 


Hence 


xWW'  2'  =  a'7WW'y 


By  denoting  the  nxn  symmetric  matrix  as 

T 

I 
t=l 

then 


2W(?W)'   =     I     w2  f(t)f'(t)   =  F 


xWW',?'    =  a'F  B-3 

Now  the  conditions  for  F  to  have  an  inverse  are 

(1)  There  are  at  least  as  many  observations,  x, 
as  degrees  of  freedom  in  the  model 

(2)  The  fitting  functions  are  linearly  independent. 

If  the  above  two  conditions  hold  true  then  by  postmultiplication  of  B-3 

by  F-1 

2    -1 
a'  -  xW  «7'F  -"- 

This  is  the  expression  for  the  vector  of  coefficients  that  minimize  the 
weighted  sum  of  squared  residuals. 
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APPENDIX  C  Correlation  and  Autocorrelation 
Given  the  tvo  variables 

h  "  *r  x 

yi  =  yl'  y2'    '    '    '    '  yn 
The  correlation  coefficient  r  may  be  expressed  as 


r  =      I   (x.   -   x)(y     -  y)    /    (n)S  S 
1=1     1  X  X  y 


where 


S^     J   (x.   -  x)2   /  n  , 


C-l 


i=l 


S2 


g  =  I  (y,  -  yr  /  n  . 

y       i=l 


If  the  variables  are  first  converted  to  standard  form  then  C-l  may  be 
expressed  more  conveniently.  That  is  if 

ui  =  Cxi  "  x)  I   Sx 
V.  =  (y.  -  y)  /  Sy 

then 

n 
r  =  I   u  v  /n  C-2 

1=1  x  1 

Tvo  properties  of  the  correlation  coefficient  are  significant 

(1)  -1  _<  r  <_  +  1 

(2)  r  =  +_  1  only  if  the  points  lie  on  a  straight  line. 

Hence,  the  correlation  coefficient  is  a  measure  of  the  strength  of  the 
linear  relationship  between  the  two  variables.   If  the  correlation  coef- 
ficient is  1  then  the  two  variables  are  directly  linearly  related.   If  it  is 
-1  then  an  inverse  linear  relationship  exists.   A  correlation  coefficient 
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of  zero  indicates  that  no  linear  relationship  exists  between  the  tvo 

variables . 

The  idea  of  correlation  betveen  two  variables  can  be  extended  to  the 

consideration  of  the  lag  correlation  of  a  single  variable.  Consider  the 

variable  x(t),  t  =  1,2,  ...,  n,  where  x(t)  represents  the  value  of  x  at 

time  t.  The  idea  of  autocorrelation  is  a  measure  of  the  correlation,  or 

the  strength  of  the  linear  relationship,  between  x(t)  and  xCt+k).   Consider 

the  expression 

T 

I       x(j)x(j-k) 

where  T  is  the  range  of  the  series.  The  expected  value  of  this  tern  is 

called  the  average  lagged  product.  Moreover,  if  the  series  has  been  adjusted 

so  that  the  expected  value,  E(x),  is  zero,  where  x  denotes  the  adjusted 

series,  then  the  average  lagged  product  is  the  autocovariance 

T 
R  (k)  =  V   x4x,  ,  /  T-l-k 
J=k+1  J  J"k 

where  R  (k)  is  the  autocovariance  evaluated  for  lag  k. 

Since  the  variance  of  a  sequence  of  numbers  which  have  been  adjusted 

to  have  an  average  value  of  zero  is  Just  the  expected  value  of  their 

"2 
squares,  E(x  ),  by  virtue  of  the  notation  Just  defined  this  can  be  expressed 

as  R  (0).  Hence,  by  C-2  the  autocorrelation  coefficient  can  be  expressed 

as 

p(k)  =  R   (k)/R   (0)  . 

XX      XX 

The  set  of  values  of  this  function  for  all  lags  k  =  +  1,  +  2,  ...  is  called 
the  autocorrelation  coefficient. 
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APPENDIX  D  Exponential  Smoothing 

The  smoothed  statistic  for  exponential  smoothing  is  defined  as 
S  (x)  =  ox(t)  +  (l-a)S^  , (x) 

o  <  I  ±  ± 

where  x(t)  is  the  observation  at  time  t 

r 

q(   is  an  undimensioned  positive  ratio  less  than  1. 
If  the  above  definition  of  exponential  smoothing  is  defined  as  single 
smoothing  then  double  smoothing  can  be  defined  as 

sJ2J(x)  =  asW(x)  +  (l-a)s|f}(x)  . 

Similiarly  multiple  smoothing  of  order  k  is  defined  as 

sW(x)  =  sf-1)(x)  ♦  (l-a)sW(x) 

The  fundamental  theorem  of  exponential  smoothing  proves  that  it  is  pos- 
sible to  estimate  the  (n+l)  coefficients  (derivatives)  in  an  nth  order  poly- 
nomial model  by  linear  combinations  of  the  first  n+l  smoothed  statistics. 

THE  FUNDAMENTAL  THEOREM  OF  EXPONENTIAL  SMOOTHING 

If  the  observations  x(t+t)  are  represented  by  the  model  (1,133) 


tkx(k  ■■'■■ 
where  t  is  the  forecast  period  then 


dt+T)  =  I  xk  x;kVk: 

k=0    t 


'(x)  =  l   (-l)k  (x<k)/k!)  ap/(P-l)!  J   jkSJ(p-l+j)!/j: 

k=0  *—n 


k=0  J=0 


where  x    is  the  k   derivative  evaluated  at  time  t. 

As  a  proof  to  the  fundamental  theorem  of  exponential  smoothing  two  vec- 
tors are  defined,  a  vector  x  which  represents  the  infinite  sequence  of  obser- 
vations x(t)  for  t  =  -,...-1,0,1,...,-.,  and  a  vector  S  with  components 


16U 


0    (t  <  0) 

oB   (t  >0)  . 

Exponential  smoothing  can  be  represented  as  a  convolution  of  the  two  vectors 
x  *  S  which  has  the  components 


(x  *  S)t  =  St(x)  =  I   x(t-j)S  =  a  I   3Jx(t-j) 


J=0      -    j=0 

Having  defined  the  convolution  relationship  for  single  smoothing,  since  the 
convolution  operation  is  associative,  multiple  smoothing  of  order  p  is  equiv- 
alent to  the  convolution  x  *  S  p  where  Sp  necessarily  has  the  components 

0  (t  <  0) 

(sW)t  = 

aV  (P-i+t):/t:(p-i)   (t  >  o) 

Therefore 

S.p(x)  =  [x(t-j)  (aV(p-l+j):)/j(p-l):  . 

1     J=o 


But  since 


c(t-j)  =  I  (-i)(x<k)/k:)  jk 
k=o    * 


the  theorem  is  proved. 

The  first  five  smoothed  statistics  can  be  written  as 
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x(k) 

s.U)  =  I  (-i)k\r-«  I   jV 
z  k=0  J=0 


'       k=0      *-     J=0 


x(k)  3 

S(3)(X)=  £(.!)*  JL-f-  [  jk(J+l)(J+2)8J 
*        k=0       K"   "*  J=0 


f,\       n      x'k'  It  » 
S.W(x)  =  I  (-l)k4r-f-  I  Jk(j+D(j+2)(J+3)6J 

1      k=o     k-    °    2=0 


■  W-  £(_i)kJL«   J  Jk(j+I)(j+2)(J+3)(J+1*)B 
k=0      k-  ^  J=0 


From  the  structure  of  these  terms  the  following  is  determined. 

(1)  The  p   order  of  smoothing  is  given  as  the  alternating  sura  of 
the  n  coefficients  in  the  Taylor  series 

xf>/k: 

(2)  The  coefficients  of  these  derivatives  are  infinite  sums  in- 
volving the  smoothing  constant. 

By  referring  to  the  closed  form  of  the  infinite  sums  (1,135). 
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Form 

Zjsj 
IjV 

v  >•  i 
ZJ6BJ 


1/1-B 


8/(1-6)' 


B(i-e)/U-s) 
b(i+hb+b2)/(i+s) 

B(1+11B+11S2+B3)/(1+B) 
B(1+26s+66s2+26b3+S  )/(l+B) 
B(l+5Te+302S2+302B3+57B  +S  )/(l+B) ' 


and  writing  the  vector  of  coefficients  as 


4° 

4" 


(t) 


xt  /n. 


the  fundamental  theorem  can  be  expressed  as 


S  =  Ma 


D-2 


vhere  M  is  an  nxp  matrix  with  elements  involving  infinite 
sums  of  cowers  of  the  smoothing  constant 


ik 


nrrrr  I* 


kDj  (i-n-.Q: 
J! 
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.  ta 
where  M.,  is  the  element  in  the  i   row  ar.d  .< 
Ik 

column  of  the  matrix  M' 
S  is  the  vector  of  the  p  smoothed  statistics  evaluate 

at  time  t.  3y  nrcessity  n*p. 
a  is  the  vector  of  the  coefficients 
With  the  fundamental  theorem  in  matrix  fcrr.  the  ti+1  simultaneous 
equations  in  the  n+1  derivatives  given  in  .1-1  can  he  solved  To-  these 
derivatives  in  terms  of  the  n-*-l  smoothed  statistics.   If  3-2  is  post- 
multiplied  by  M~  the  result  is 

a  =  Sjf1. 
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APPENDIX  E 

The  purpose  of  this  appendix  is  to  present  the  computer  programs. 
The  programs  are  referred  to  by  figure  numbers  in  the  following  manner. 
E-l     The  first  phase  of  the  program  for  general  exponential 

smoothing  -  "INCONT" 
K-2     A  model  for  a  linear  trend  plus  a  superimposed  sinusoid 

(included  in  Phase  I) 
E-3     Phase  II  -  "RAY" 
K-U     Phase  III  -  "MATINV" 

E-5     The  subprogram  to  calculate  the  h  vector  -  "HVEC" 
E-6     Phase  IV  -  "FORCST" 
E-7     Phase  V  -  "PLOTTER" 
E-8     The  subprogram  to  plot  the  observations,  forecasts  and 

errors  -  "PLOTS" 
E-9     Monitor  cards  used  for  phasing  in  the  PR-155  system 
E-10    The  program  for  calculation  the  autocorrelation  function 
E-ll    The  subprogram  to  remove  the  trend  from  the  time  series  - 

"TREND" 
E-12    The  program  to  calculate  the  power  spectrum 
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414  READ! 1 .929) (X( I ) ) 
Ptah  IN  THF  PFRICD  NTOU 

RrAD(  1.1)  (NTC'J  ) 
6?^  IFtNNN-?  )  77,88,88 
88  STO° 
■77  GO  TC(  31 ,66)  ,NNN 


BF  TRIFP 


FIG. 


F-l 
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»3.14*8/TCU) 
. 14*B/TCU) 


'3    CONTINUE 
v  N  =  4 
N  =  4 
C  I  (  1  >  a  ?  2  0  7  a 

rr<2)»-l?.? 

r  I  <  3 1 = o .  o 

CI(4)=C.O 

DC3I=1 ,NTCU 

TCU*NTCU 

8  =  1 

C3  =  X( I  )*SIN( 2. 

C  r ( 3 ) =C I ( 3!+C3 

C4  =  X( I )*CCS<2.*3. 

3  c:  (4i  =  ci  (  4i+o 

C I ( 3 )  =  (2. /TCU)*CI t  3  ) 

C!  (4)  =  (2.  /TCUI*CI  C») 

ad,1  1=1. 

A (?,] )=0. 

6 (3,1 )=C. 

A  (  A  ,  1  )=1. 

CHK (11=0.0 

CHK (2 1=1.0 

CHK (31=1.0 

CHK  (4)  =0.0 
RFAD     IN    THF     TRANSITION    MATRIX 

TM( 1,1 )=1.0 

Tv( 1,21=0.0 

TNM 1,3 1=0.0 

T"'(  1,41=0.0 

T  W  <  2 , 1  1  =  1.0 

T*M  ?,?)  =  ]  .0 

T  M  (  ?  ,  3  )  =  0  . 0 

TM(2.4)=0.C 

Tf(3tl l=O.C 

Tf(3i21"0.0 

Ty(3,3)=CCS(2.*3.]4/TC'J> 

TM<3»41=SIN(2.*3.14/TCU) 

TV  (  4,  i  )  =  0  .0 

TV (4, 2 1=0.0 

TW(A,3)=-<;IN(?.»3.14/TCU) 

TM(4,4 )=C05<  ?.*3.]4/TCU) 

O^TORi 8 
f-f,    CCNTINIIF 
8  I  8    v:R  I TF (  7  1  NO ,NTCU» N •  (  A  (  I  •  1 1  ,  I  =  1  » N  1 

•  {XU)»I-l»ND)»(CI<l)»I-l»N)»NNNi 
1KN»  (CHK(  I  1  ,  I  =  1,N1 

DOB  171  =  1  ,N 
817    V  R  I  T  F  (  7  )  (  T  v  (  I  ,  J  )  ,  J  =  1  ,  N  ) 

NNM=NNN+1 

GO    TO    625 

rND 


FIG. 
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MOM**      FXFO  FORTRAN. . ,0§». . .PAY 

J  NTFGFRRNP 

PI  MFN'C  I  ON  A  t   (  1  r  .  1  ) 

n  t  v  r  v  c  i  « v  r  *•  (  o  ,  9  j  .  r  y  l  (9,0) 

p  r vcvc ; ONA ( 9 . 1  !  .  D ! 1  ,  9  ;  ,C]  (  9.Q  )  , F(< 

D I M  E  NS I ON  C  0 1  (  9  ) 

DIMENSI0NX(5UG) .CI (9) .A1I9.1) »H<9.1] 


1 
111 
11? 
990 
900 
91  - 

6A? 
'l  A  I, 


?.l) .TM<9»9! .AF<9»1 ) 
A) 


STAPTING 

9.7 


DIMENS10NA 
FORMAT  (  14  ) 
FORMAT  (  n  X  ,  ]C'F1  1 
rOPMAT  ( 1X.I4) 
rORMAT( IX, 6! 10) 
"CPMATM X.F3C.7) 

■CP«AT(irx .121 
,  r;on«T  (  r-F  1  =.  .  5  ) 
6?<j  FOPMAT  (  1  X  ,  ]  ?F°.^  ) 
75  F0PMAT(3(  X.F2  >.A] 
2?  FORMAT (5F 20. 4) 
2 A  FORMAT! 1UX »5F20. 4 ) 
624  FORMAT! 1X.18H  F  INVERSE 
REWIND7 

pfw  r  \ir>5 
RFWIND6 

MPMsl 
D  MD  =  ^ 
'"*♦*#-}<  A      rift.'     MOOFL       IS      PFAP      T  A' 

C»##**RETA     Ic     !N!TT2L!ZFn    tq 

Aft      A  9  F  T  A  -  ;    .  7  (j 

"=""  A  '3    (  7  IND.NTOU  »N»  !  A(  I  •  1  )  .  I  = 
:  .  <X(  I  )  •  1  =  1  ..ND)  .  (CI  t  I  )  ,  I*1,N) 
IMN. (CHK(  I  ) ,I»1 .N) 
D08 171=1 ,N 
817    RFAO    (7) (TM( I ,J) ,J«1 ,N) 
T  M  =  N 

A  p  =  ( 1  .  /  T  N  ) 
r       av'"    THF     INITIAL    VALI 
Pp9^6 1=1 . M 

9  * '.  ■",T(n='"i(r) 

r»*nncnvc  THF  IMITIAL  VALUFS 
O*hhh<0F  THF  FITTING  FUNCTIONS 
rT77I  =  7  ,N 
77  A0(I»1)=AII»1) 
( »##**STAPTING  HERE  TUT 
C*»***»FCAt  CUt  ATFD  WITH 
67  9"FTA= ( AHETA ) **AP 
C*****PfINITIALIZE  THF  VFCTOR 
C**#**CF  FITTING  FUNCTIONS 
P04?  1  =  1  .N 
f*I  (  I  I  ^r^.'   (  I  1 
h 7     AA(J.1)=&0(T.11 


HFRF 


>N  ) 
NNNi 


F    THF    CONSTANTS 


VATr> ix     I ? 
NEW    !?FTA 


F  I  r, . 
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■id 


Di'Tf    THF    f    i/JTRIX    F"° 
"^^ri       f.:ir>     ^rtr     pcta 

[  r  1  ,  N 

J  s  1    .  « 

!  .J! =0.0 

.  j )  = :; .  i. 

A  » 1  .  0 


0 1  =  1  .  N 

, I  )=A(  I  ,1  ) 
I  =  1  ,  M 


"1  ( 
TN 

'=1 


c*»* 

c»** 
c»»* 

ft 
A 


I.J)=A ( !»1 )*B( 1 »J )*(ncTA ) 

TINUE 

31=1 »N 

3J=^»N 

(  I  ,J)=C]  (  I  ,J)*:JETA 
K-l  )66  ,16,66 

C<    THr    CONVERGENCE    CF    THE 
ATRIX    AT    EACH    5CTH    ITERATIt 
S    ROUTINE    CHFC*:S    UNTIL     IT 
r,c    aw    ELEMFNT    WHICH    I r- 
c  t  A  R 1 1  1 2  F  n 

CNT-50.  )  I  6. A/.  ,1  ft 

1  =  1  tN 


C»»* 
c**» 

^  *  s  * 

c*#« 
c»*« 


TPI 


Lj  =  r- 

M*F< 
M=C1 
10=' 
RATI 
TIN'I 

'■'  A  ! 

V      A  | 

X    »E 


J) 

.J! 

M/DFNM 

.10000] >80»80.16 


IAS 

)Cg 

.  =  1 

DCH 
:K  ( 

■  ( I 

.'Q  T 

.  (A 
Tft 


"  TH 
PFTF 
TIT 

PFI 
71=2 

-! 

7J  =  1 
I  ..J) 
•  J )  = 
TC(ft 

1(1, 

r  i  =  l 

TT(ft 


[T    FROM    THIS    L""n    OCCURS 
T(jp    ||  rvrt'T'.    OF     THF    F    N»A- 
THF    CONVERGENT    CRITERIA. 

OCCURS     THr     F     "AT^IX     I  s 
r>Y     ^YMJiFTRY    ANH    THI-     PPCGRAY 
SEE    ~F    tljf    RANGE    Cr     F1ETA 

cc 


:vrRF: 


^c" ( J,  I ) 

F  (  J  ,  I  ) 

!  n viv  .  M  .n">  ,  (     CI  (  I  I  •  !  =1 

l  )  .  I  =  l  .  M  ,  (  X  ,  !  )  .  I  =  1  , «. 


.N)  • 

M  ,  i  A  i 


)  (  T"-'(  I  . 


:  1  ,  M 


(CCNTINt  FD) 
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'■">  I  T  F  (  r-,  )  V  ,  (  ,1  C  (  I  ,  I  )  .  I  =  1  .  N  )  .  r  «  c  T  A 

i>f>68!s]  .N 

r. r  l»,<»ITF{MfFU»J)»J«1»N) 

00481=]  ,N 

/.  »  WITFt^l  (FK(  I»J)  ,J«1  ,N) 

C#**»*PARAMETRIC  0*i  PETA 

IF(ARETA-0.70)102»45>45 

45  RNR=RNR+] 

I  F(  OMR-tvQM  )46,i  4,  t  /4 

]  4  STOP 

J  f\  2      A  n  r  T  A .  =  A  ^ETA+0.20 

C«»»**RCPLFNISH  THP  VFCTCR  Cp 

c*#***fitting  functions  and  use  a  mfw 

C#*#**R  PTA 

PC  T"  ft  7 

C-»»#**THF  PROGRAM  PRANCHFS  TC  HERE  WHEN 

C***»«TH?  CONVERGENCE  CRITFRIA  IS  NOT  "FT 

C*****AMD  NORMAL  ITERATION  CONTINUES. 

16  P05U  =  1  iN 

fW  =  J 

D050J=M,N 

FK  (  I  ,J)rFK(  I.JI+CIK!  .J) 

5^  F(  !  ,J)=F  (  I  .J)+C1  (  I  .J) 

C*****UPDATF  THC  VPTTCR  "F  FITTING  FUM- 

C*«*««CT!CN<  USING   THF  TRANSITION  ,va- 

C##***TPIX 

r>03]  1  =  1  »N 

?]  AF( I , ] 1=0.0 

00901=]  .N 

D090J=]  .N 

9"  AF( I,1)=TK( I.J)*AA( J,l )+AFt 1,11 

rr;f:ri  =  '  ,n 

DC1C0J=1 ,N 

AA(I,1 1=AF(I,1) 

100  MI,]  )=AF(  1,11 

C»#*#*THIS  ROUTINF  CONVERTS  THE  VICTOR 

r>**»*"F  FITTING  FUNCTIONS  EVALUATFD  FOP 

{■#♦**»«  fiych  TIME  ORIGIN  TO  THF  VcrTPP 

C#»**»P\/ALUATFn  rCP  A  MOVING  TIT  ORIGIN 

[>063I  =  ] «N 

rv  =r\-tv   ;  I  ) 

I  F  (  C-  I  .  )  6?  .  55  .ft' 

55  A( I  ,1 ) =  -A (  I  , ]  ) 

63  CONTINUE 

C  SAVF  F( 1 1  1  OR  FCRbCAST 

!F(K-1 162,??, 62 

??    D062I=! ,N 

A](I,])rAA(I,l) 

62  ""'!^f 

FIG.   r-?  (CONTINUFD) 

17>t 


BFTA*CSETA* 

CNT«CNT+J . 

<  «K+1 

G  °    T  "    99 

STOP 

Ff>'D 


:TA 


-1     (CONTINUED) 
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7?  7 
4  44 

777 

8P8 
8P9 


MSN'S 
I  IT 

[5  I  >.-' 
CI'-' 

r>  TV 
■  P  I  M 

n  i « 

F^r? 


66 

67 
33 


1THF 
FOR 
FC" 
PNR 
N  P  N 
RFW 
R  F'.»> 
pes 

rice 

Q  FA 

DC3 

RFA 

DC6 

D06 

Ct  I 

NM1 

DO 

ALL 


ERR 
SIC 

MO 
SI" 
s  I" 
SIC 
M/ 
TIF 
TI9 
T  (  / 
CEF 
T(/ 
T(/ 


FX 

N" 
NFK 

NV? 
VV1 

M  A  C 

NC( 

f  n 

F  !  3 
//! 
FIC 
//l 
//2 


CO    FORTRAN. 


vat  i,\y 


.10) 
.10] 


7  )  » H  (  9  »  1  ) 
0  ]  »V(9)  iF(<>,9  ) 

1AH     TIJC      (-)     VC<"TC 


.3  ) 

0X»39H  THE  VARIANCE  : 

I  E  N  T  S  I  S       ! 

X  ,  1 1"  F  I NVFR5E  .    ) 

0X»6H  r;FTA=.F9.6) 


P7 

P  5 


44 
11  = 
!  ) 
99 
1  =  1 
'11 
55 
.J 
»N 
DO  66 
6669  C.  (N.J 


144  ^ 


*  5  5  0 

9990 


VI  Jl 

V  (  N  ) 
DO  ( 
tPJ< 

CIP] 
DC 
C(  I 
C(  I 


:  ) N»  (  AO (  I  • 1  )  .  1  =  1  .N)  .O^FTA 

= !  »  N 

5) (F( I.J) »J  =  ] -N ) 

=  1  ,N 

5)<FK(I»Jt»J=l»N) 

l.N 

1  .N 

) =  F (  I.J) 

-1 

6>.  L=1.N 

!1»1) 

4;   J=]  .N^I 

r< i  .J+1 )/all 
1 ./all 

9  i   I  =  1 , NM ] 

+  1 

=  C(  I P  1  » 1 ) 

5-  J=1.NM1 

)=C( IP1.J+1)  -  CIP11*V( J ) 

)=(-!. 0)*CIP1 1*V(N) 

6'..  J  =  1.N 

) =V( J) 


FIG. 


•r-h 
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r*  n  6 1 « 1  •  \i 

\)C16J=\  »N 

16 

F(I  «J)='C(  I.J) 

C  CALCULATE  THE  H  VECTOR 

CALL  HVfcC ( F.AC.N.H) 

C»*#*#WRI  T! 

dF.TA 

WRITFI 3.8891CBETA 

C  *  *  *  ->H 

;v.'P  I  TF  F  I  NVERSF. 
WR1 TF ( 3 . aRP ! 
DC3?3!=1 .N 

■m 

'■"?!  T^l  3  .'.44  )  (  F  (  !  .  J)  . 

J  =  l .Nl 

C    WRITF  CUT  thf  H  VFCTCR 

W  R  I  T  r  (  *  ,  7  7  7  ) 

DC62fi!=l »N 

WRITFI3.789) (H(I»1) ) 

676 

WRITF(2.789)  (h(  I  .1  )  ) 
DC7I=1 »N 

7 

WRITFI 7) <F< I.J). J=li 

N) 

WRITE(7) (H( I.1)»I=1. 

N) 

C  »  *  *  *  '■ 

tFCRN  THE  V  MATRIX 
r>  ?  1 4 1  =  ]  » N 

n«i  4j  =  l  ,v 

14 

V?( I »J)=C.C 

r)cni=i  »n 

DC11K=1  ."1 
001 1J»] .N 

1  1 

V?(K. I )=c (K»J)*FK( J, 
DC13I*]  .N 
DC1 3J=1 »N 

I )+V?( K.I ) 

1  3 

Vll I.Jl'O.C 

DO] 7  1  -  '  »N 
DC1 2K»1  ,N 

""1  7J-1  .N 

i  7 

VI (K.I )=V?(K.J)*F( J, 

WRITFI ^.777) 

!  )+V  (K.I) 

177 

WRI TF( 3.789) V1 (K»K) 

n  |.IPsPMP  +  ] 

If (RNR-NRN166.68.68 

6R 

STCP 
F  ND 

FIG. 


F-4     (CCNTINl.'r  D) 
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c,|nip"i  it  I  mp    HVFC(F»AC»N 
"I'Tv-ci  2.N r  (  Q  .  9  )  •  A j  (  q  .  ] 
7PO    i-  f^'.MT  (  r  1  ■» .  5  ) 

"-  ?     H  (  I  ,  ]  )  s  0  .  L 

DCS^  J  =1  »  N 
DC543J=1  » ,\ 
HF( I»l >=F( I »J)*AC(J»1 1 

54'     u  (  I  ,1  )=H(  I  ,1  )+HF(  i  ,1  ) 

s?n    rCNTINI'F 

•JfTCPf,. 


.lJ(  c 


jyFr 

.!Jt;  [9,1 


FIC 


178 


rypn  r  o.  P  T  p  A  V  , 


,  .FORCST 


222 

1  "7 

1^8 

1  ^1 

1 

1  OR 
0  7? 

67? 

105 
1(4 
106 


r  ^t 
r.JM 

P!V 
15 1 M 

1  .TV, 
D  I  !>' 
DIM 
Frp 

1  Tur 


C***# 


FOR 
F 
1 
FOR 
C0R 

FOR 
o  r-.i 

PC 

or 
VPM 

P  r^ 

i  (  x  ( 
WPJ 

W  I 

*wr  i 

*VAl 

WRI 

0^7 


FUPV 
C  IQIVj 

SIGN 

SIGN 

9)  . 

ION 

51  CM 

(  // 

RrC 

(  1  X 

T(  IX 

T  (  1 X 

=  .!~4 

T  {  IX 

T(72 

RFCA 

T(F2 

T(l< 

T(  18 

T!24 

-it 

17 


F(9 
CCI 
X<5 

CIF 

Al  ( 

TTM 

/]X 

AST 

,F2 

.15 

.11 

.?) 

,?C 

H 

ST 


,F1 

H  S 


.9) 
(9) 
00).H(9»1) »CI(9) . 

(  9  )  •  C  i  L  (  9  ) 

9,1  ) . A ( 9  ,  ]  ) 

(9.9) 

,31H  THF  VAPIANTK  OF 

S«,F1 5.6/1  HI) 

0.7.F20.4] 

H  NFW  FORFCAST    ) 

H  "CHrL  NC  .  =  •!?•  "OX  «6H 

H***THE  F^ROR  IS***        ) 
PERIOD         OBSERVATION 

ERROR     CUM  ERROR  ) 
) 

.3) 

5  .  '■  •  F  1  5  .  2  .  F  I  b  .  2  . r  1  5  .  2  ) 
UM  OF  SQUARE  OF  FRROR= >F15. 3  J 


r   ( 
i ) . 

TF( 

TE( 

Tr 

UF 
TF( 
1  =  3 
5(6 


6  )  N)  **'  M 
1 1 1  ,  N 

5  )  M  D  . 
2,1  OP 
OUT  T 
OF  F>- 
3  »  1  0  ] 


, N . nn  .  (  C I  I  I ) » I  =  1  »N 1 . ( A 1  (  I  .  1  ) » I  =  1 »N ) ! 
[) )  ,ASFTA  ,MN 
(X( I)»I=1.ND) 

) 

HF  MCDEI  NUMBER  A;*'p 

TA  FOR  THIS  F  MATRIX 

JMN.AFETA 


■"A  T  : 


OTA 

p  r  ^ 
sr "  fi  vr 
« VAT 


D(7 

D(7 

F  T 
"IX 


1  ( 
,  M 

)  C7!  I 
I  (H(  I 

HF  To 


( T»J)»J«1 .M) 

i  J  )  »  J  *  1  .  N  1 
,  1  )  .  I  a  I  .  N  ) 
ANSPOSF  O1-"  TMF  TRANSITION 


FIR. 
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Sti"X=f«C 

i?o    tt*>'(  I  ,J  )=TM(  J,  I) 
r»**-**>'A<r    rHE    FORECAST 

',•0  1  tf  (  3  ."i ! 

Z=0.f! 

n  o  1 1  i  I  = 1  » v 
ZM=CI( I)*A1( I tt) 
711     Z*Z+ZN 

VRITFI51Z 

FOeX(f )-Z 

r»**»«CALC!ILATc    THF    fUM'l 


ATIVr 


51 11 
Slil 


IMA+CO 


C***#*r ALCHL  ATE    Tac    sum    OF    T'-J,~    On^FRVATIONS 

SUMX=SUMX+X< K) 
C    UPDATE    THF    CONSTANTS 

D0449I  =  1  iN 
449   CIL( I )=0.0 

D094M =1»N 

f  IF(  I)«TTM( I f J)*CI ( J) 

1/.4     Ml  (  I)=r  IL  (  T  l+CIF(  I  ) 

D09?fi  T  =  "^  »N 
o?A    <"  f  {  7  )  =r  1  (_(  I  I  +H(  T  »1  )  *  CP 

VR  IT'r(3»104)  (K»X(K I  >Z  »FR.SU"A ) 

WRITR(?!  104)  (K.XK)  »7  >FR) 
193     CONTINUE 

i'.'RITFI  3.198) 

WRITFI3.106) (SUM) 

'.•'PI  TF(2»107)  SUM, BETA 
C*****CALniLATE  THE  VARIANCE  OT  TUT  FORECASTS 

VAP  =  <-.!"'/c."vX 

WITF  [3>???)VAR 

l?VOsPMR  +  ) 

J  p  (  [JWp-MPM  )  ftf.  .  r,P.  ,',P 
ftfl     r,  TOP 

I'Mp 


FK 


(  CON  I  INU,:D 
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VCN'SJ  FXEO.    FORTRAN.  » »09»  •  t  tPLCTTEfi 

T  NTFGFRRVP 

DIKFNSIONX(500) . FCST ( 5  00) 
RFWINr>5 
is  fv'  f? »  r 

"l)M=1 

66    PFAP(^)ND,(X(  I )  .  !*1.ND) 

ps    FCST(K)=Z 

CALL    D  L  C  T  5 ( X , F  C  S T  ) 

R  N'R  bRM  R+  ] 

IF(RNR-NRN)6&»68»68 
6fi    STOP 


PIT-.      F-T 
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.^V(r»               rXr"    FrjOTPiN,  ,  ,0o,  .  ,  ,plC 

■ 

c.t  itr?"l  IT  !  MTD[    "Tr,|Y.7l 

riT'/rvr  J;  VY(  1  )  ,7.  I ]  1  •vD  (  1  ??) 

7 
i 

F0PVAT(1X»132A1  ) 

r  "ii"f.T  1^1  ,?  T'> ) 

4 

FT"  "'AT  I1X.1VI1H+]  ) 

i 
5 

FC9.VAT  (  1 1'1  ! 

FCRMAT!//ll-X«40H            TI«r    SERIES. 

F  OP  FrAc>T  ,  FRRCR  ) 

6 

FCR"AT I//1CX.13HX=TIME    SER!rP»10X, 

1  r.M"  =  c"D-rA?T»     1  CX  .7H^-r::'r'nr}  ) 

•7 

11 

F  "DM  AT(  //T"X,AH"AX  =  .  I  A.=X.AHVIN=, 

'  ■>  ,^,A|-f/'Lr=  .  !  •*  ! 

WP I TF( 3,3) 

QCADI1  .1  )WS]  ,(k'S?,MS3,WA,MR,MAXt»'I'ii»NPTS 

WOITFI 3,4 ) 
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In  many  industrial  and  economic  situations  a  series  of  discrete  ob- 
servations are  taken  on  a  uniform  time  scale.  This  sequence  of  obser- 
vations is  called  a  time  series.  The  object  of  this  investigation  is  to 
develop  a  method  for  forecasting  these  time  series. 

If  the  series  of  observations  is  not  merely  random  it  can  be  described 
as  being  the  sum  of  two  components.  One  component  is  the  process  vhich 
generates  the  time  series  and  the  other  is  the  variation,  or  superimposed 


